МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ РАДИАЦИОННО-КОНДУКТИВНОГО ТЕПЛООБМЕНА В ПЛОСКОМ СЛОЕ ПОГЛОЩАЮЩЕЙ И РАССЕИВАЮЩЕЙ СРЕДЫ Текст научной статьи по специальности «Физика»

УДК: 536.24

Математическое моделирование радиационно-кондуктивного теплообмена в плоском слое поглощающей и рассеивающей среды

В.С.Виноградов , О.Н.Третьякова , Д.В Хакимов .

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

В ряде технологических процессов возникает необходимость использовать сверхтонкие пленки, нанесенные на разогреваемые поверхности [1]. Требуется определить поле температур и распределение интенсивности излучения или тепловых потоков в системе "пленка-подложка". Это приводит к необходимости решения задачи радиационно-кондуктивного теплообмена (РКТ). Задача РКТ в общей постановке представляет собой связную систему интегро-дифференциальных уравнений переноса излучения и теплопроводности, решение которой сопряжено с большими математическими трудностями.

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

Основными методами решения уравнений переноса излучения, записанных в дифференциальной форме, являются: метод сферических гармоник ( рм ), метод средних потоков (СП метод), метод дискретных ординат Чандрасекара ( SN ), метод собственных функций Кейса, метод экспоненциальной аппроксимации ядра и некоторые другие методы [2-5]. Рассмотрим отдельные методы с целью выяснения возможности применения для решения поставленной задачи. 1) Метод сферических гармоник [2]. Запишем уравнение переноса излучения в следующем

виде

1 —I —I 3 ^ i

--^ + = ~kvIv +avÈv +3\ïv (x, ¡, n')Iv (x, ¡, t)d^\

° ni riir ¿Lit j

vv

ïv(x,¡,л ) = j/v(x, ¡)d0 , где ¡л = ¡¡' + -V1 - ¡2 V1 - л2 cos(0 -0), в -азимутальный угол,

0

л = cos0 - направляющий косинус угла, Iv - спектральная интенсивность излучения, <^v -

коэффициент поглощения, pv - коэффициент рассеяния, kv = а, + - коэффициент ослабления,

В, - функция Планка, У, - индикатриса рассеяния.

Решение ищется в форме ряда по полиномам Лежандра

^ +1

I, (х, t) = х —— л, (х, t) рп (м),

4+

+1

(х, г) = 2+ | К (х, М г)Рп ЫФ .

-1

Умножим уравнение переноса на Рт (м) и проинтегрируем по [—1,+1]. Индикатрису рассеяния тоже разложим по полиномам Лежандра

у, (х и) = X 24 +1 Уп,(х) рп (¿0, 4+

Уп

-1-1

(х) = 2+|уу (X,м) рп .

—1

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

да л . 1

у, (х, и, м ) = X —Г" Уп, (X)Рп (М)Рп (и').

п=0 2

Интегральный член тогда преобразуется в

Р V,......ч, ..' Р ^ 2п +1

4 (х м, и У,(х и, t=р- X~4—У",(х)**т (х't)Р" (м).

4+ —1 2 п=0 4+

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

+ 2^+7 + , \ Рт С"Ж

1 + т д*Щ — ХУ т + 1 д*т _

с, дг 2т +1 дх 2т +1 дх тУ тУ У У —1 т

1 +1 где kmV = kv — 1 ЯУтг , ^ = ^ =ау , Уо, (х) = 2, ^ (х) = Г Р0С")Рт = 0 , при т Ф 0

2 —1

В случае Р1 приближения (базис {Р0(и), Р1(и)} = {1, и}), применительно к стационарному полю излучения, получаем систему уравнений относительно плотностей потоков объемного падающего Л*, и результирующего излучений. Совместно с граничными условиями,

определяется приближенное значение интенсивности излучения

I,(х,и) = [*0,(х) + 3и* 1у(х)] = 4+ Л(ь) + 3иЕ,(Ь)] = Л(Ь) — и^Ь!]. 4+ 4+ 4+ ат„

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

2)При исследовании радиационного и комбинированного теплообмена наиболее широкие возможности составляет метод средних потоков [2] и его модификации. Основная идея метода средних потоков заключается в замене интегро-дифференциального уравнения переноса энергии излучения, записанного в виде

л у

МТ = -1V + (1 -«V Я + « 11V (7, М)У, (Н, М)йм системой из двух дифференциальных уравнений, получающихся последовательным

интегрированием этого уравнения по М е [—1,+1] с весом 1 и М , соответственно. = ^/п + п -

спектральное альбедо однократного рассеяния.

Используя представление о встречных потоках, приходящих из полусфер с противоположными знаками О = ±2л в сечении тv получаем

-Т [ Е+ + Е— ] + (1 — )(— + —— ) = 4(1 — ,

ат„

V

-7 [П+ + П - ] + (1 — «у (1 — ))(Е+ + Е— ) = 0,

йТ;

где Л*; - плотность потока падающего излучения, Еу - плотность потока результирующего излучения, П у - диагональная составляющая тензора излучения, ¿>у - функция формы индикатрисы рассеяния.

Вводим коэффициенты

— — П+ П

,„,,+ IV „„— IV 1+ V ?— V

ту = — • т„ = — • I„ =-• / =-

V + 5 V т т — 5 V + 5 V — •

Е Е — —

V V IV IV

Получаем систему двух дифференциальных уравнений относительно — и —— с коэффициентами m— , т—, 1+, I—

й - [—7 + ] + (1 — «V )(— (7) + П— (7)) = 4(1 — «V)К (Тv):

dТv т+ (Тv ) т— (Тv )

-Т [С (Т — (7)+1— (Т —— (т )]+(1—«V (1—¿V М^+тН+—тЧ ] = 0.

dТv ^(Х) т— (Тv)

Используя конечное число членов разложения уУ (р, р ) в ряд по полиномам Лежандра, по заданному начальному распределению температур находим поле значений интенсивностей излучения. Затем вычисляются т+ , т~, IV , 1~, ЛУ+, ЛУ , . Найденное значение плотности

потока объемного падающего излучения Л*У = ЛУ + ЛУ, подставленное в уравнение сохранения энергии, позволяет находить температурное поле, используемое для вычисления в последующих итерациях.

3)3адачи РКТ в различных постановках с использованием различных модификаций метода матричной прогонки решались авторами работ [4-7]. В работах [6,7] были рассмотрены задачи РКТ с использованием точных уравнений переноса излучения при произвольном характере индикатрисы рассеяния.

Для решения практических задач в качестве первого приближения можно использовать упрощенную модель решения задачи РКТ в плоском слое нерассеивающей среды, предполагая, что вклад рассеянной энергии мал (область температур 300-500 К). Рассмотрим следующую задачу. Имеется плоский однородный слой вещества (подложка - 1), на которую напылили пленку (граница - 2). Геометрия задачи представлена на рис.1.

Рис. 1. Геометрия задачи РКТ Для моделирования нашей задачи рассматривается плоский слой серого излучающего и поглощающего вещества с диффузно отражающими и излучающими границами. Граничная поверхность 1 поддерживается при постоянной температуре Т1, на границе 2 задан радиационный тепловой поток. Границы имеют степени черноты и ¿>2 . Соответствующие уравнения и граничные условия примут вид [3]

д1 (г, р) п 2аТ4 п

Р-г-+ 1(г, Р) =-, (1)

дг ж

n 2aT4 1

I(0, p) = s1 П ° 1 + 2р? f I(0,-p)p dp л i

I (T0,-p) = £.

n 2оТ4

л

1

■2p2d J I (To, p)pdp ,

(3)

Н

где Р = сое в - направление между распространением излучения и осью х, г0 = |а( x)dx -

оптическая толщина, «(х) = а1,0 < х <Н -И0,а(х) = а2,Н -Ь0 < х <Н - коэффициент поглощения,

Н - толщина подложки с напылением, - толщина напыления, р1, ^ - коэффициенты диффузного отражения левой и правой границы, п - средний показатель преломления среды. Формальным решением является (+ прямое направление, - обратное направление)

-т/ T 1 п 2сТ 4(т') -(T-T'V I + (т,p) = I + (0,p)e /^JV (T )-

Р

л

']e

dT

Г(т,p) = I"(т0,p)e /"-J—/^т

,(г0-г% п 2сТ V) -(г-г) 0Р

Заменим в (5) р на - р

л

I- (r,-p) = I- (0,-p)e

(т0^ т° 1 п2аТ4(г')п -(т-т)

о p

л

]e

'dr

Подставим в (4) т = т0, а в (6) т = 0

- У (• I

I+ (T0,p) = I+ (0,p)e /p+j-p[

1 rn2стТ4 (т )

p

- У г 1

I - (0,-p) = I - (0,-p)e /p + f-[

0 p

л

0 1 n2аТ4 (т )

]e

_(т> -т )/

1 dz

]e7 pdт'.

л

Подставляя (7) и (8) в граничные условия (2) и (3), получаем

I + (0) = a1 + b1a2 1 - bb '

I - (т ) a2 + Ъ2a1 I (т0) =

1 - bxb2 Ъ = 2р^Ез(г0), b2 = 2р^Ез(т0),

п 2аТ4

ax =

л

2р1 Л

1 т» 2_т4

п аТ (т )

0 0

л

e /pdidp

2 ti4

n аТ2

a2 = + л

2pd JJ

1 т0 „2 т-г4

n 2аТ (т )

_(т -т )/

0 0

л

! dт dp

(5)

(7)

(14)

(4)

(6)

(8)

(9)

(10)

(11) (12)

Г —У

где Е3(т0) = I не Мйм - интегро-экспоненциальная функция.

0

Для определения температуры решается следующая краевая задача дТ (х, t) Я 5 2Т (х, t) 1 дqrпd (Т (х, t), х) 1

dt cp р dx2 Г(x,0) = T0 , T (0, t) = T1,

cp р

dx

cp P

■f.

Я

dT dx

= arad\

4 x= H

где f - функция, задающая внутренние источники тепла (постоянная).

(16)

(17)

(18)

(19)

Алгоритм решения:

1) Решаем уравнение теплопроводности без учета излучения для получения начального распределения температур Т[x,t,0] .

2) Определяем I+(Т[х,¿,0],х,н),I— (Т[х,¿,0],х,н) .

+1

3) Определяем qrad (T[x, n,0], x) = I (T[x, t ,0], x, ju)dju .

cp P

dx

4) Решаем уравнение теплопроводности

dT[x, t,1] = Я d2T[x, t,1] ^ 1 dqrad(T[x, t,0],x) ^

dt cp р dx2

T [ x,0,1] = T0, T [0, t ,1] = T1,

cp P

-f.

Я

dT [ x, t ,1]

dx

= a (T [ x, t,0], h ).

5) По найденному распределению температур определяем новое поле интенсивностей I+ (T[x, t,1],x, u), I— (T[x, t,1],x, u) .

6) Повторяем итерационный процесс до нужной точности.

Численные расчеты были произведены с помощью программы Maple 7.Коэффициент поглощения моделировался двумя функциями: кусочно-линейной и экспоненциальной с большим провисанием. Данные о коэффициенте поглощения взяты по экспериментальным данным из книги [8].

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

1лверхнее Лверхнее * /- л \

1 г 4пк(Я)

серое л л

верхнее н^нее Янижнее

Ja(2)d2

2 — 2 ^

вернее нижнее Яижнее

Я

d2

(21)

x=H

x=H

Я — Я

верхнее нижнее

± ^ (Як — Я,—1),

к=1 Як

(20)

где R - количество разбиений отрезка длин волн, которые учитываются в расчете поглощения. Вычисления производились для видов материалов пленки: алюминий и никель. На рисунках 2,3,4 представлен вид натурального коэффициента поглощения , по

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

Рис. 2. Вид натурального коэффициента поглощения для алюминия,

= 1095.58

Рис.3. Вид натурального коэффициента поглощения для никеля,

(Хс

= 126.22

Рис.4. Вид натурального коэффициента поглощения для

1

подложки из кварцевого стекла,а

серое

= 8.19

Коэффициенты поглощения подложки а1 и пленки а2 существенно различны а1 << а2.

Для моделирования коэффициента поглощения как функции координаты мы использовали кусочно-линейную аппроксимацию а(х) = а150 < х < Н -Ь0,а(х) = а2,Н -Ь0 < х < Н и

/ \ % /1 а 1ч х.

экспоненциальную аппроксимацию а(х) =а1 + х ехр(1п(-%—)—) • Z - параметр вида

Н Н

функции. В численных расчетах непрерывность экспоненциальной функции сказалась на скорости работы программы. Разложение кусочно-линейной функции в ряд Фурье для ускорения процесса вычисления не принесло успеха, так как малое количество членов ряда (< 100-200) дают физически не правильные (отрицательные) значения для коэффициента поглощения. Для случая кусочно-постоянного коэффициента поглощения расчет выполнялся примерно 40 минут. Экспоненциальное же а снижает время расчета до 8-10 минут (вычислительная машина Репйит-4 2.53 ГГц, 512 Мб ОЗУ). На рис. 5 представлен вид кусочно-линейной и экспоненциальной зависимости коэффициента поглощения для толщины напыления ^0 ~ 31-25мкм (при кол-ве разбиений М=320) для алюминия.

Рис. 5. Вид кусочно-линейной (1) и экспоненциальной (2) ^=100) функций поглощения. Температурные поля и выходящий радиационный поток. Для дальнейших расчетов поля температур в исследуемой пленке и подложке данные были взяты из справочников [8-10]. е1 = 0.85, 82 = 0.85, р? = 0.08 , Т0 = 300К , Тх = 400К , ср = 500Дж /кг • К,

Х = 0.74 Вт / м • К, р = 2210 кг/м3 , п = 1.5, / = 107.

Время расчета t = 10 секунд. Толщина пластины с пленкой Н = 1 см. М=320 - количество разбиений отрезка по толщине [0,Н]. На рис.7.-16 результаты расчетов для q = 0,когда решается только уравнение теплопроводности без учета излучения, представлены красным цветом, для q =1, 2, 3 , когда решается связанная задача РКТ с использованием одной, двух, трех итераций соответственно, представлены синим, зеленым и черным цветом.

Начальное распределение температур ^=0 с.) задано, как показано на рис.6.

Рис.6 Начальное распределение температур. Для алюминия с использованием кусочно-линейной аппроксимации коэффициента поглощения

а(х) := <

8.19 х <

319 32000

1095.58

319 , при оптической толщине слоя г0 = 012 и коэффициенте

32000

$ х

диффузного отражения р2 = 0.95 [11] получены распределения температуры по толщине слоя для разных моментов времени, что представлено на рис.7-9.

Рис. 7. Температурное поле для алюминия при t =3 с.

Рис. 8. Температурное поле для алюминия при t = 7с.

Рис. 9. Температурное поле для алюминия при t = 10 с.

Видно что прогрев в слое не равномерный. Температурное поле с учетом лучистого переноса тепла в теле, характерно возрастает на пленке (при t=10 с., Т2 - Т % 50 ). Анализ графиков показывает, что итерационный процесс уже к 4-5 итерации выходит на приближенное решение уравнений РКТ. Сходимость явно зависит от количества разбиений временного интервала, а так же от шага времени.

Рис. 10. Зависимость от времени выходящего результирующего радиационного потока с границы пленка-вакуум для алюминия при х = Н.

Рис. 11. Динамика нагрева правого неизотермического края х=Н для алюминия. Для никеля использовалась функции поглощения экспоненциального вида

„ х 100 (46531.30570)

а(х ) := 8.19 + х е , Z=100. Оптическая толщина г = 0.84,

коэффициент диффузного отражения р2, = 0.66 [11].

Температурная динамика для никеля показана на рис.12-14.

Рис. 12. Температурное поле для никеля при t =3 с.

Рис. 13. Температурное поле для никеля при t =7 с.

Рис. 14. Температурное поле для никеля при t =10 с.

Рис. 15. Зависимость от времени выходящего результирующего радиационного потока с границы пленка-вакуум для никеля при х = Н

л т2(я)

¿20-

400-

380-

360

Ш"

320-

300

Рис. 16. Динамика нагрева правого неизотермического края x=H для никеля.

Выводы. Даже с принятыми упрощениями видна сложность задачи нахождения решений уравнений РКТ. Выполненный расчет показал, что нанесением тонкой металлической пленки можно влиять на распределение температур в теле и на его радиационные характеристики в целом. Из сравнения рис. 10 и 15, 11 и 16 , соответственно, видно, что для никеля выходящий радиационный поток и температура правого края ниже, чем для алюминия. Основной задачей в каждом конкретном случае является выбор материала пленки с нужными коэффициентами поглощения С и отражения р . Влияет так же и толщина пленки. Динамика прогрева, зависимость лучистого переноса тепла от коэффициента поглощения, отражения и толщины пленки представляют актуальную задачу, которую авторам предстоит в дальнейшем решать.

[1] Резник С.В. Математические модели радиационно-кондуктивного теплообмена в материалах тепловой защиты многоразовых транспортных космических систем. // Инженерно-физический журнал. - 2000, Том 73, №1.- с.11-25.

[2] Рубцов Н.А. Теплообмен излучением в сплошных средах.- М.: Мир,- 1984.- с.74-77, 85-88.

[3] Оцисик М.Н. Сложный теплообмен.- М.: Мир,- 1976.- с.295-302.

[4] Формалев В.Ф.Тепломассоперенос в анизотропных телах. Обзор.//Теплофизика высоких температур. - 2001,т.39, №5.- с.810-832.

[5].Третьякова О.Н Численные методы математического моделирования задач сложного теплообмена. Деп в ВИНИТИ № 6310-В89.-185 с.

Литература

[6].Третьякова О.Н. Анализ термонапряженного состояния пластин из технического стекла. Дисс. к.ф.-м.н., М.:1983.-194с.

[7].Липовцев Ю.В, Третьякова О.Н. Конечно-разностное решение одномерной нестационарной задачи радиационно-кондуктивного теплообмена.// Инженерно-физический журнал.-1986, т.51,№ 5.- с.840-847.

[8] Золотарев В.М., Морозов В.Н., Смирнова Е.В., Оптические постоянные природных и технических сред. Справочник. - Л.: Химия, 1984.

[9] Кошкин Н.И. и Ширкевич М.Г. Справочник по элементарной физике. М.:1972,- с.37,74, 85-86.

[10] Прядко А. Журнал "625" Светотехнические материалы. 3/2005. - http://rus.625-net.ru/625/2005/01/prjadko.htm

CВЕДЕНИЯ ОБ АВТОРАХ

Виноградов Владимир Сергеевич, доцент кафедры физики Московского авиационного института (государственного технического университета); к. ф.-м. н.; Телефон: 158-48-24, e-mail: vin 42@rambler. ru

Третьякова Ольга Николаевна, профессор кафедры физики Московского авиационного института (государственного технического университета), к. ф.-м. н.; Телефон: 158-86-98, 8(916)542-03-18,

e-mail: tretiyakova_olga@mail.ru

Хакимов Дмитрий Викторович, студент факультета "Прикладная математика и физика" Московского авиационного института (государственного технического университета); Телефон: 793-37-65

e-mail: 7933765@mail.ru