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

УДК 536.3:535.34

К РАСЧЕТУ ЛУЧИСТЫХ ТЕПЛОВЫХ ПОТОКОВ В ПРЯМОУГОЛЬНЫХ ОБЛАСТЯХ МЕТОДОМ ДИСКРЕТНЫХ ОРДИНАТ

А.В. САДЫКОВ

Нижнекамский химико-технологический институт (филиал) ФГБОУ ВО «Казанский национальный исследовательский технологический университет»

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

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

Радиационный теплообмен является основным в большинстве технологических трубчатых печей нефтехимической промышленности. При заданном поле температуры, радиационных свойств излучающей среды и ограничивающих поверхностей расчет лучистого теплообмена в стационарном случае сводится к решению интегро-дифференциального уравнения переноса излучения (УПИ) при соответствующих граничных условиях. Для решения УПИ применяется множество различных методов. В задачах лучистого теплообмена широко используется метод дискретных ординат (МДО), предложенный Чандрасекаром [1]. Метод получил развитие в работах как зарубежных, так и отечественных ученых. Преимуществом этого метода является то, что он удачно объединяется с алгоритмами, основанными на применении контрольных объемов. В частности, метод был использован автором статьи в ряде работ [2-4 и др.].

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

У 0,6

Стенка 2

0,0 Стенка 1

Рис. 1. Система координат

а,0

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

вид

д/ (г, О) д1 (г, О) Ц-1-+£- '

дх

ду

а1ь (Т) - (а + Р)1(г, О) + [ I(г, О')у(О', П)ЛП'. (1)

4 %

© А.В. Садыков

Проблемы энергетики, 2016, № 3-4

Здесь ц, £ - направляющие косинусы; I(г, О) - интенсивность излучения,

зависящая от положения и направления; а, в - коэффициенты поглощения и рассеяния

среды соответственно; у(О ', О) - индикатриса рассеяния; Д (Т) - интенсивность

излучения черного тела при температуре Т.

Граничное условие к уравнению (1), с учетом диффузного излучения и отражения, имеет вид

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

В МДО УПИ (1) заменяется системой дифференциальных уравнений относительно интенсивности излучения вдоль ограниченного количества выделенных

направлений {т = 1, 2, ..., N0}. Эти направления задаются угловыми координатами {ц„, £т ; т = 1, 2, ..., N0}, равными величине проекции единичного вектора направления на оси координат 0х и 0у соответственно. Исходя из числа

выделенных направлений, различают ^-приближение (N0=4), ^-приближение (N0=12) и другие. Таким образом, получаем систему дифференциальных уравнений относительно интенсивности излучения Iт вдоль каждого из этих направлений:

где - угловые весовые коэффициенты. Индексы т', т обозначают направления

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

соответствующим направлению $ телесным углом ДО.

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

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

Показателем асимметрии является коэффициент g\, который заключен в диапазоне - 1 <1. Для изотропного рассеяния = 0 (фт 'т =1).

Граничное условие (2), характеризующее излучение раскаленных стенок и отражение лучей, падающих со всех направлений, в МДО для различных стенок аппроксимируется следующим образом:

(2)

для таких направлений О , для которых п ■ О > 0 . Здесь О ' - направление

у(0) = 1 + £1С08 0 .

ф - =1+ЯГи Ц -+£ '+Л Л '1.

N

Im = sIb (Tw ) + - Z wm' Vm' I7m' % T, m =1

при x=0 (стенка 3 рис.1) |m > 0 и Vm' < 0; при х=а (стенка 4 ) |m< 0 и Vm' > 0;

N0

Im = sIb (Tw) + - Z wm' \%m' \I,

(5)

'm'=1

при у=0 (стенка 1) Е,т > 0 и Е,т- < 0; при у=Ь (стенка 2) Е,т < 0 и Е,т- > 0.

В выражениях (4)-(5) вторым членом в правых частях описывается отраженный поток лучистой энергии, при этом суммирование ведется только по направлениям падения лучей.

Дискретные аналоги уравнения (3) получаются интегрированием его по контрольному объему вокруг точки (/,у), показанному на рис. 2. В результате получаем систему алгебраических уравнений:

(1т(/+1, у ) — 1т(/—1, у)) + (1т(/, у+1) —1т(/, у-1)) _ у — V;, у + у ;

рг,у = а/,у1Ъ (Т/,у )*/,у ; V/,у = (а/,у + Р/,у К',у ;

S>,J =

j Gi, j

N

4%

Z wm 'фmmIm (i,j) ; G/,j 4BiAj ;

m' =1

Aj = 0,5(yj+1 - yj); Bi = 0,5(хг+1 - Xi).

>

(6)

j+1 j

j-1

1

i

i

~r

"T"

p.................О................Ö

6.................6................<>

77777777777777777777777777777777777?'.

1 г—1

i

i+1

Рис. 2. Разностная сетка

Решение системы уравнений (6) для каждого из ординатных направлений находится методом итераций. В каждой итерации используется метод покоординатной

прогонки. Допустим, что искомая величина I ^ в узловой точке (/, у ) связана со

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

1т(/, у) = ш1т(/, у +1) + (1 — 1т(/, у—1) = ш1т(/+1, у) + (1 — ш)1т(/—1, у), где ш — интерполяционный коэффициент.

Формулы для численного решения алгебраических уравнений получаются следующим образом. Рассмотрим направление интегрирования, для которого цт>0 и ^т>0 и индексы возрастают. Значения интенсивностей 1т(/1 у—1) и 1т(/1_\ у) будем считать известными, которые назовем опорными значениями. Из соотношения (7) выражаем

интенсивности 1т^ , ,) через опорные значения и подставляем в уравнение

системы (6). В результате приходим к следующей прогоночной формуле:

_ ^тА]1т(1-1,]) +^тВ11т((,]-1) + ю№,] + ] )

Пусть теперь 0, ^ < 0. Опорными значениями интенсивностей считаем

1т(,-1,]) и 1т(1,,+1) . Из соотношения (7) выражаем 1т(1,,-1) , Iт(,+1,,) и поДставляем в

уравнение системы (6). В результате получаем формулу для значения интенсивности в средней узловой точке I

т(',])

га

т(ьЛ

НтА]1т(г-1,]) ^тВ11т(1,] +1) + га№',] + ])

А га „ „ пщ , + ЦтА, ---^тВ1

^ -1 1 - га

Для других значений угловых координат и ^ получаются аналогичные выражения.

Итерационный процесс для решения системы (6) совместно с дискретными аналогами граничных условий (4), (5) реализуется по следующей схеме.

1. Задается начальное приближение для I ,) в узловых точках области.

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

3. Для каждого ординатного направления (цт, (т=1, 2, 3, ..., N0) обходятся все

узлы разностной сетки и вычисляются значения интенсивностей ^ по

соответствующей прогоночной формуле с учетом знаков и ^ .

4. Полученные значения интенсивностей I^ принимаются за начальные

значения для следующего шага итерации и происходит переход к пункту 2. Итерационный процесс продолжается до выполнения условия

тах

]

п ..и+1 %]

<8,

где п - номер предыдущей итерации; 5 - заданная малая величина;

Nо

фг- ■ = ^ ^„¿тц ]) - объемная плотность энергии излучения.

т=1

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

N0 N0

^х ^ ^т^т^т ' ^у ^ 1 ^т^т^т .

т=1 т=1

В работе [5] приведены результаты численного решения задачи лучистого теплообмена в двумерной прямоугольной области, заполненной изотропно-рассеивающей, поглощающей и излучающей средой (рис. 3), зональным методом, ©Проблемы энергетики, 2016, № 3-4

методом сферических гармоник (в Р3-приближении) и точное решение. В указанной работе рассмотрены три случая: 1) объем с изотропно рассеивающей средой с черными границами; 2) объем с изотропно рассеивающей средой и с серыми границами; 3) объем с поглощающей средой с черными границами. Здесь приведены результаты численного решения этих задач в 52-, 54-, ^-приближениях МДО с помощью разработанных нами программ. При этом использовались конечно-разностные сетки 10^20, 20^20. На всех рисунках приведены результаты, полученные с использованием сетки 20x20. Численные исследования показали, что дальнейшее дробление сетки в случае однородной среды не приводит к уточнению результатов. Направляющие косинусы и весовые множители такие же, что и в работе [5].

Рис. 3. Область интегрирования. Индексы параметров соответствуют номерам стенок

Итерации прекращались при достижении максимального отличия значений объемной плотности лучистой энергии в двух последовательных итерациях не более 0,1%. Результаты расчетов приведены в безразмерном виде. Значения поверхностных плотностей радиационных тепловых потоков отнесены к интегральной энергетической светимости абсолютно черного тела соответствующей

температуры ^ = оТ4 = ^ / ^ ). В качестве линейного масштаба принята длина

стороны квадрата.

Изотропно рассеивающая однородная среда

Предположим, что излучающей стенкой является 1 -я стенка = 1), для других границ приняты условия =~4=0. В расчетах критерий Шустера 5с=1,

оптическая толщина среды г=1.

На рис. 4 показаны графики изменения безразмерных поверхностных плотностей потока результирующего излучения на стенке 1 от ее края к середине для значений ее степени черноты: ^ = 1; 0,5; 0,1. Результаты, полученные с использованием 54-, 56-приближений хорошо согласуются с результатами расчета по зональному методу. Для абсолютно черных границ (8 =1) отличия не превышают 2% для 56-приближения, 7% - для 52-приближения.

О 0,1 0,2 0,3 0,4 ^

Рис.4. Распределения безразмерных поверхностных плотностей лучистых поттоков к ст. 1: -----Р3 [6],----- В4;-----®, I, А - зональный метод [6]

Для серых границ (8 =0,5; 8 =0,1) для В4-, ^-приближений наблюдается хорошее соответствие с зональным методом, при этом кривые для этих приближений практически совпадают. Поэтому на этом рисунке для этих случаев приведена кривая только для В4-приближения.

При уменьшении степени черноты стенки точность В2-приближения ухудшается. При 8 =0,5 отличия достигают почти 15%. Однако следует отметить, что в топочных устройствах степень черноты обычно меняется в пределах от 0,6 до 0,9.

Ухудшение точности В2-приближения при уменьшении 8 объясняется тем, что увеличивается доля отраженного излучения, и в граничных условиях преобладающим становится слагаемое, учитывающее падающее на стенку излучение, а собственное излучение стенки становится незначительным. Ввиду того, что в МДО угловое распределение интенсивности излучения аппроксимируется конечным числом угловых интервалов, то анизотропии рассеяния и отражения в низших приближениях этого метода не учитываются.

Графики изменения относительной объемной плотности энергии излучения и в зависимости от безразмерной координаты у для сечений X =0,5; 0,3; 0,1 при т=1 показаны на рис. 5. При этом границы считались абсолютно черными.

й 2,0

1,5 1,0 0,5 0

и

2,0 1,5 1,0 0,5 0

Результаты расчетов, полученные с использованием 52-, 54-, ^-приближений сравниваются с результатами, полученными зональным методом и в Р3-приближении метода сферических гармоник (МСГ) [6]. Использование 52-приближения приводит к завышенным значениям вблизи горячей поверхности и к заниженным - вблизи холодной стенки. Как видно из рис. 5, погрешности являются наибольшими вблизи средней линии (X =0,5) и убывают при приближении к боковым стенкам (X =0,3; 0,1). Следует также отметить, что вблизи боковой стенки (X =0,1) 54-, 5^приближения дают почти одинаковые результаты.

Однородная среда без рассеяния излучения

Здесь даются сопоставления результатов расчетов с аналогичными результатами, полученными при следующих идеализированных случаях: 1) границы двумерной прямоугольной области полностью поглощают излучение среды, но сами не излучают (~ = ~ = ~ = ~ = 0> £1=£2=8з=84 = 1 - абсолютно черные холодные стенки); 2) в среде отсутствуют центры рассеивания (5с=0); 3) мощность излучения среды ©Проблемы энергетики, 2016, № 3-4

поддерживается равной единице (¿~=1). Оптическая толщина среды принимает значения: т=0,1; 1; 10.

Распределения безразмерных поверхностных плотностей потока результирующего излучения на стенке 1 при трех оптических условиях показаны на рис. 6.

Ч

0,7

0,5 0,3

0 0,2 0,4 0,6 0,8 X б

1,2

1,0

0,8

0,6

Г = 10

, — _--- •V

\ ч

/ _ - . — _ . ч

/ ^— 1 \

' / \\

// \\

// \\

0 0,2 0,4 0,6 0,8 X

Рис.6. Распределения безразмерных плотностей результирующих потоков излучения: --точное решение [5]; — - — ■ — ■ - - ,-----

Решение, полученное в 52-приближении, дает значения q , значительно отличающиеся от точного решения во всех точках для оптических условий г=0,1; г=10. Б4-, Б6 -приближения дают значения с более высокой точностью, причем наибольшие отклонения наблюдаются в центре поверхности. Результаты расчетов, полученные с использованием этих приближений, почти совпадают. Например, при т=0,1 максимальные отличия результатов не превышают 0,001. При таком выборе масштаба кривые практически совпадают. Поэтому на рисунке кривая для ^-приближения не приведена.

Во всех рассмотренных случаях 56-приближение обеспечивает наилучшее согласие с точным решением, при этом отличия результатов расчетов от точного решения не превышают 6%. В проведенных численных исследованиях интерполяционный коэффициент га =0,5.

Преимуществом МДО, как уже было отмечено выше, является сочетаемость с алгоритмами, основанными на применении контрольных объемов. Кроме этого, важным преимуществом МДО является то, что для перехода от низшего приближения к более высокому достаточно изменить в расчетной программе значение N0, массивы ©Проблемы энергетики, 2016, № 3-4

Ч

ím, Wm (m=1, 2, 3, ..., N0). Опыт использования МДО и МСГ при решении задач

сложного теплообмена показал, что МДО экономичен как в плане затрат машинного времени, так и в плане необходимого объема оперативной памяти компьютера.

^-приближение непригодно для расчетов излучения в чисто поглощающей среде. Причиной расхождений между численными результатами и точным решением является «лучевой эффект». Такая погрешность обусловлена конечной дискретизацией угловой переменной в уравнении переноса излучения. Излучение может распространяться только вдоль фиксированных дискретных направлений. Излучение от изолированного источника может не восприниматься некоторой точкой, если эта точка и источник не лежат на каком-либо из дискретных направлений. В задачах с чисто поглощающей средой и малым числом ординатных направлений точки области слабо связаны с соседними точками. Члены, описывающие рассеяние внутрь, связывают данную точку с ее соседними точками, в результате чего уменьшается погрешность метода. «Лучевой эффект» можно ослабить, увеличивая число ординатных направлений. В большинстве задач, имеющих практическое значение, рассеяние нужно учитывать. При уменьшении степени черноты стенки точность ^-приближения ухудшается. S2 -приближение рекомендуется использовать при значениях 8 в пределах от 0,6 до 1.

Литература

1. Чандрасекар С. Перенос лучистой энергии. М.: ИЛ, 1953. 431с.

2. Садыков А.В., Вафин Д.Б., Садыкова Д.А. Зависимость тепловых и аэродинамических характеристик трубчатых печей от расположения ярусов веерных горелок // Известия вузов. Проблемы энергетики. 2014. №11-12. С.3-10.

3.Садыков А.В. Методика расчета сопряженного теплообмена в трубчатой печи производства водорода в рамках дифференциального метода // Известия вузов. Проблемы энергетики. 2013. №7-8. С.3-11.

4. Садыков А.В. Влияние подсоса воздуха на сопряженный теплообмен в трубчатой печи производства водорода // Известия вузов. Проблемы энергетики. 2014. №1-2. С.37-44.

5. Fiveland W. A. Discrete - ordinate solutions of the radiative transport equation for rectangular enclosures // Trans. ASME: J. Heat Transfer. 1984. V. 106, №4. P. 699 - 706.

6. Ratzel A., Howell J. Two - Dimensional Radiation in Absorbing - Emitting - Scattering Media Using the P N Approximation // ASME Paper No. 82 HT 19, 1982.

Поступила в редакцию 15 марта 2016 г.

Садыков Айдар Вагизович - канд. техн. наук, доцент, декан факультета управления и автоматизации Нижнекамского химико-технологического института (филиал) (НХТИ), г. Нижнекамск. Тел.: 8(8555)392314; 8(917)8624162. E-mail: sadykov_av@mail.ru; sadykov@land.ru.