УДК 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: [email protected]; [email protected].