CC BY
9ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ / PHYSICAL & MATHEMATICAL SCIENCES
УДК 536.3:535.34
К РАСЧЕТУ ЛУЧИСТЫХ ТЕПЛОВЫХ ПОТОКОВ В ТРЕХМЕРНОЙ ПРЯМОУГОЛЬНОЙ ОБЛАСТИ МЕТОДОМ ДИСКРЕТНЫХ ОРДИНАТ
TO CALCULATION OF RADIATION HEAT FLUXES IN THREE-DIMENSIONAL RECTANGULAR ENCLOSURES BY METHOD OF DISCRETE ORDINATES
©Садыков А. В.,
канд. техн. наук, Казанский национальный исследовательский технологический университет, г. Нижнекамск, Россия, sadykov_av@mail.ru
©Sadykov A., Ph.D.,
Kazan National Research Technological University, Nizhnekamsk, Russia, sadykov_av@mail.ru
Аннотация. В работе метод дискретных ординат применяется для решения задачи переноса излучения в трехмерной постановке.
Приведены результаты, полученные с использованием S8 — приближения, проведено их сравнение с данными других авторов.
Abstract. In work the method of discrete ordinates is applied to the solution of a problem of transfer of radiation in three-dimensional statement.
The results received with use of S8 — an approximation are given, comparison with data of other authors is carried out them.
Ключевые слова: интенсивность излучения, радиационный теплообмен, контрольный объем.
Keywords: intensity of radiation, radiation heat exchange, control volume.
В камерах радиации технологических трубчатых печей тепловые потоки формируются в результате совместного действия всех механизмов переноса энергии: лучистого, конвективного и кондуктивного. При температурах, характерных для трубчатых печей нефтехимической промышленности, лучистый теплообмен является основным. Расчет лучистого теплообмена сводится к решению интегро-дифференциального уравнения переноса излучения (УПИ) при известном поле температур. В задачах лучистого теплообмена часто используется метод дискретных ординат (МДО).
Первоначальный вариант МДО как развитие двухпотокового приближения Шустера и Шварцшильда для плоской геометрии был предложен Виком [1] ив дальнейшем развит Чандрасекаром [2]. Основное предположение метода: представление интегрального члена в УПИ по квадратурной формуле. МДО получил развитие в работах зарубежных и отечественных авторов. Он удачно объединяется с алгоритмами, основанными на
научный журнал (scientific journal) Т. 4. №5. 2018 г.
http://www.bulletennauki. com
применении контрольных объемов. Метод был использован автором статьи в ряде работ [34].
Рассмотрим УПИ для прямоугольной области, показанной на Рисунке 1. Для излучающей, поглощающей и рассеивающей серой среды уравнение имеет вид
+ а1ЛТ) - (а + р)/(г,П)+А \ЦТД')уф,ТШ . (1)
дх ду & 4тг (4->л)
Здесь ц, п — направляющие косинусы; I(г, О) — интенсивность излучения,
зависящая от положения и направления; а, р — коэффициенты поглощения и рассеяния
среды соответственно; у(О', О) — индикатриса рассеяния; (Т) — интенсивность излучения
черного тела при температуре Т .
Граничное условие к уравнению (1) с учетом диффузного излучения и отражения имеет
вид:
I(r, Q) = s/j (T) + - Цn Q'[I(r, n')dQ' (2)
Л nQ'<0
для таких направлений О , для которых п • О > 0 . Здесь О' — направление
падающего излучения, О — направление испускаемого излучения; в — степень черноты граничной поверхности; р — отражательная способность поверхности; п — единичный вектор внутренней нормали к границе.
В МДО угловое распределение интенсивности излучения аппроксимируется значениями вдоль ограниченного количества выделенных направлений 8т {т = 1, 2, ..., N0}
в каждой точке пространства. Для трехмерного поля излучения эти направления задаются набором угловых координат {цт, ^т, Пт; т = 1, 2, ..., N0}, равных величине проекции единичного вектора направления 8т на оси координат 0х, 0у, 0z соответственно. Количество
направлений в Sn — приближении равно ^=и(и+2), где п — порядок приближения. Интегро-дифференциальное уравнение (1) приближенно заменяется системой дифференциальных уравнений относительно интенсивности излучения вдоль каждого из этих направлений:
dI dl dl ß No
+ + = aIb - (a + ß) Im + Г" 2 mmIm , (3)
dx dy dz 4л m '=i
m=1,2,3,...,M
Vo
где — угловые весовые коэффициенты. Индексы т', т обозначают направления
падающего и испускаемого излучения соответственно. Угловые весовые коэффициенты Wm связывают между собой интенсивности излучения вдоль различных направлений. Эти коэффициенты численно равны площади на поверхности единичной сферы, отсекаемой соответствующим направлению 8т телесным углом АО. Коэффициенты и направляющие косинусы удовлетворяют условиям:
X ^ = ; ^2т +^2т +Л2т = 1.
т=1
Индикатриса рассеяния, характеризующая рассеяние лучистой энергии мельчайшими частицами сажи во всех направлениях записывается в упрощенном виде:
у(0) = 1 + g cos0 .
Коэффициент фт,т , учитывающий анизотропию рассеяния при квадратурном представлении интегрального члена, определяется по выражению
ф - =1+gГ и и -+£, '+л л -I.
т m т o1Lr~ 1т 1т J
Показателем асимметрии является коэффициент gi, который заключен в диапазоне — 1 <gi <1. Для изотропного рассеяния gi = 0 ( фт,т =1).
Граничное условие (2), характеризующее излучение стенок и отражение лучей, падающих со всех направлений, для различных стенок аппроксимируется следующим образом:
Р ^
1т = ^b (Tw ) + Wm' U т'К' (4)
Л т =1
при х = 0 для ит < 0 и Um > 0, при х = Lx для ит > 0 и и™< 0;
Р No
I =£I (T ) + РУw ,\t М - (5)
т b \ w J / * m|^mpm V /
Л т '=1
при у = 0 для ^т' < 0 и ^т > 0, при у = Ly для ^т ' > 0 и ^т < 0;
Р N
1т = (Tw ) + ^' |Лт ' К' (6)
% т' =1
при z = 0 для Пт' < 0 и Пт > 0, при Z = Lz для Пт' > 0 и Пт < 0.
В выражениях (4)-(6) вторым членом в правых частях описывается отраженный поток лучистой энергии, при этом суммирование ведется только по направлениям падения лучей.
Область интегрирования покроем конечно-разностной сеткой. Уравнение (3) проинтегрируем по контрольному объему, показанному на Рисунке 2. В результате получаем следующее конечно-разностное уравнение относительно интенсивностей излучения в узловых точках:
итАр (1Ьт - 1т ) + ^тВр (£ - ^ ) + ЛтСр (£ - С ) = Fp-у pIPm + Sp . (7)
Здесь Ар, Bp, Cp, ур — функции координат и оптических свойств среды в точке P расчетной области D; Fp — источниковый член в точке P, зависящий от интенсивности собственного излучения среды; Sp — источниковый член, определяющий вклад в
интенсивность излучения в данном направлении вследствие рассеяния по другим направлениям.
научный журнал (scientific journal) Т. 4. №5. 2018 г.
http://www.bulletennauki. com
Рисунок 1. Система координат. Рисунок 2. Контрольный объем.
Значения угловых координат и весовых множителей для первого квадранта приведены в Таблице.
Решение системы уравнений (7) совместно с граничными условиями находится итерационным методом. В каждой итерации используется метод покоординатной прогонки. Для получения дополнительных уравнений предполагается, что интенсивность излучения в центральной узловой точке контрольного объема связана со значениями в соседних узловых точках следующими выражениями:
1Р = ъ1ат + (1 - ъ)1ст = ш^ + (1 - ш) 1ьт = ъ1ет + (1 - ш) 1т, (8)
где ю — интерполяционный коэффициент. При ю=0,5 получаем «ромбовидную» схему, а при ю=1 — «шаговую» схему или схему «против потока». Формулы для решения алгебраических уравнений получаются следующим образом. Рассмотрим направление интегрирования, для которого Цт > 0, ^т > 0, Цт > 0. Из уравнения (7) с помощью соотношений (8) исключим интенсивности в узловых точках Ь, с, / В результате приходим к следующей прогоночной формуле:
Jp ^т^р^т + ^>тВр1т +ЛтСр1т + ш(-^р + ) р + МгаАр + ^>тВр +ЛтСр
Значения интенсивности в узловых точках Ь, с, f затем вычисляются с помощью соотношений (8). Для других значений угловых координат Цт, ^т, Цт получаются аналогичные выражения. Например, если выполняется расчет в отрицательном —
направлении, то соответствующее рекуррентное соотношение для Im получается заменой на — и Iem на I^ . Аналогичная замена делается для отрицательного п — направления. Если
оба направления отрицательны, то делаются обе замены и т. д. Эти замены являются результатом использования уравнения (8) в уравнении (7) при нахождении интенсивности излучения.
Таблица.
УГЛОВЫЕ КООРДИНАТЫ И ВЕСОВЫЕ МНОЖИТЕЛИ ДЛЯ ПЕРВОГО КВАДРАНТА [51
m Ординаты Веса
Цт ^т пт Wm
S2 — приближение
1 0,5773503 0,5773503 0,5773503 1,5707963
S4 — приближение
1 2 3 0,2958759 0,9082483 0,2958759 0,2958759 0,2958759 0,9082483 0,9082483 0,2958759 0,2958759 0,5235987 0,5235987 0,5235987
S6 — приближение
1 2 3 4 5 6 0,1838670 0,6950514 0,9656013 0,1838670 0,6950514 0,1838670 0,1838670 0,1838670 0,1838670 0,6950514 0,6950514 0,9656013 0,9656013 0,6950514 0,1838670 0,6950514 0,1838670 0,1838670 0,1609517 0,3626469 0,1609517 0,3626469 0,3626469 0,1609517
S8 — приближение
1 2 3 4 5 6 7 8 9 10 0,1422555 0,5773503 0,8040087 0,9795543 0,1422555 0,5773503 0,8040087 0,1422555 0,5773503 0,1422555 0,1422555 0,1422555 0,1422555 0,1422555 0,5773503 0,5773503 0,5773503 0,8040087 0,8040087 0,9795543 0,9795543 0,8040087 0,5773503 0,1422555 0,8040087 0,5773503 0,1422555 0,5773503 0,1422555 0,1422555 0,1712359 0,0992284 0,0992284 0,1712359 0,0992284 0,4617179 0,0992284 0,0992284 0,0992284 0,1712359
Итерационный процесс для решения системы уравнений (7) совместно с дискретными аналогами граничных условий (4), (5), (6) реализуется по следующей схеме:
1) Задается начальное приближение для интенсивности излучения в узловых точках расчетной области D;
2) В граничных узловых точках области производится расчет интенсивностей излучения по выражениям (4), (5), (6). В первой итерации границы считаются черными, а слагаемые, описывающие рассеяние излучения во внутренние ячейки, равными нулю. В последующих итерациях используются заданные значения степени черноты стенок и учитываются члены, описывающие рассеяние излучения на частицах сажи;
3) Для каждого ординатного направления (ц т, ^т, Пт
) (т=1, 2, 3, ..., N0) обходятся все узлы разностной сетки и вычисляются значения интенсивностей по соответствующей прогоночной формуле с учетом знаков ц т, ^т, Пт ;
научный журнал (scientific journal) Т. 4. №5. 2018 г.
http://www.bulletennauki. com
4) Полученные значения интенсивностей Iтр принимаются за начальные значения для
следующего шага итерации и происходит переход к пункту 2. Итерационный процесс продолжается до выполнения условия:
max
peD
Ф'П -Фp
n+l
Ф p
<5.
N
где п — номер предыдущей итерации; 5 — заданная малая величина; фр = ^ ~^>т1тР —
т=1
объемная плотность энергии излучения в точке Р.
После завершения итерационного процесса компоненты вектора плотности результирующего потока излучения в направлениях х, у, ъ в узловых точках расчетной области определяются суммированием по всем направлениям:
No No No
Чх = X ^mWmI!p i Чу = X ^mWmI!p i 4z = X ЛЛ^ .
m=1
m=l
m=l
УПИ решено для единичной кубической области (Хх =1 м, Ly=1 м, Х2=1 м), показанной на Рисунке 1, с использованием S2, S4, S6, S8 — приближений МДО и проведено сравнение с численным решением работы [6], полученным зонным методом. Ограничивающие область стенки считались черными 8=1. Излучаемая мощность была единичной при х=1, у=1, ¿=0 и нулевой на остальных стенках. Среда находится в равновесии с излучением. Коэффициент поглощения среды а = 1,0 м-1, коэффициент рассеяния Р=0 м-1.
Угловые координаты и весовые множители такие же, что в работе [5]. Сумма весов равна 4п, значения ординат симметричны. В расчетах использовалась «ромбовидная» разностная схема второго порядка точности. Численные решения получены на сетках 5 х 5 х 5 и 11 х 11 х 11.
Результаты численного решения, полученные в S8 — приближении, приведены на Рисунках 3, 4, 5. На Рисунках 3 и 4 показана излучаемая мощность на плоскостях ¿=0,5 м и ¿=0,9 м в сравнении с решением [6]. Плотности лучистых тепловых потоков на поверхность при ¿=1,0 м приведены на Рисунке 5. Максимальные ошибки наблюдаются при ¿=1,0 м в углу при х=0 и х=1.
научный журнал (scientific journal) Т. 4. №5. 2018 г.
http://www.bulletennauki. com
Рисунок 3. Излучаемая мощность на плоскости z=0,5 м.
научный журнал (scientific journal) Т. 4. №5. 2018 г.
http://www.bulletennauki. com
0,2 0,4 0,6 0,8 x, м 1,0
Рисунок 4. Излучаемая мощность на плоскости z=0,9 м.
научный журнал (scientific journal) Т. 4. №5. 2018 г.
http://www.bulletennauki. com
Рисунок 5. Результирующие потоки тепла на стенку z=1,0 м.
Преимуществом МДО, как выше было отмечено, является сочетаемость с алгоритмами, основанными на применении контрольных объемов. Кроме этого важным преимуществом МДО является то, что для перехода от низшего приближения к более высокому приближению
достаточно изменить в расчетной программе массивы p,m, ^т, Пт, wm (m=l, 2, 3, ..., No). Применяемые для расчета лучистого теплообмена в топках печей метод сферических гармоник (МСГ) [7, 8] и зональный метод [9] более трудоемки. Кроме того, в нерассеивающих средах МСГ дает некорректные результаты [7].
Таким образом, исследования показывают, что Ss — приближение МДО дает приемлемые результаты для прямоугольной области. Достигнуто согласование результатов расчетов с численными экспериментами других авторов.
Список литературы:
1. Wick G. C. Über ebene Diffusionsprobleme // Zeitschrift für Physik. 1943. V.121. №11-12. P. 702-718.
2. Чандрасекар С. Перенос лучистой энергии. М.: ИЛ, 1953. 431 с.
3. Садыков А. В. Методика расчета сопряженного теплообмена в трубчатой печи производства водорода в рамках дифференциального метода // Известия вузов. Проблемы энергетики. 2013. №7-8. С. 3-11.
4. Садыков А. В., Вафин Д. Б. Неравномерности обогрева реакционных труб и распределений температуры продуктов сгорания по глубине технологической трубчатой печи // Тепловые процессы в технике. 2014. №8. С. 349-355.
5. Lathrop K., Carlson B. Discrete-Ordinates Angular Quadrature of the Neutron Transport Equation. Technical Report. LA-3186. Los Alamos: Los Alamos Scientific Laboratory, 1964.
6. Larsen M. E. Exchange Factor Method and Zonal Formulation for Analysis of Radiating Enclosures Containing Participating Media. Ph.D. Thesis. Austin: University of Texas, 1983.
7. Шигапов А. Б. Радиационный теплообмен в топках котлов // Известия вузов. Проблемы энергетики. 2014. №1-2. С. 27-36.
8. Шигапов А. Б., Гирфанов А. А., Ширманов М. В. Расчет радиационных тепловых потоков к стенкам топки котла в P5 - приближении метода сферических гармоник // Вестник Казанского государственного технического университета им. А. Н. Туполева. 2012. №1. С. 1823.
9. Кулешов О. Ю., Седелкин В. М. Математическое моделирование локального результирующего теплообмена в экранированных топках // Тепловые процессы в технике. 2012. Т. 4. №3. С. 118-124.
References:
1. Wick, G. C. (1943). Über ebene Diffusionsprobleme. Zeitschrift für Physik, 121(11-12), 702-718.
2. Chandrasekar, S. (1953). Perenos luchistoi energii. Moscow, IL, 431. (in Russian)
3. Sadykov, A. V. (2013). Method for calculating the conjugate heat transfer in a tube furnace for hydrogen production within the framework of the differential method. Izvestiya vuzov. Problemy energetiki, (7-8), 3-11.
4. Sadykov, A. V., & Wafin, D. B. (2014). Uneven heating of reaction pipes and temperature distributions of combustion products along the depth of the technological tube furnace. Teplovye protsessy v tekhnike, (8), 349-355.
5. Lathrop, K. D., & Carlson, B. G. (1964). Discrete ordinates angular quadrature of the neutron transport equation. Technical Report. LA-3186. Los Alamos, Los Alamos Scientific Lab.
6. Larsen, M. E. (1983). Exchange factor method: an alternative zonal formulation for analysis of radiating enclosures containing participating media. Austin, University of Texas.
7. Shigapov, A. B. (2014). Radiation heat exchange in boiler furnaces. Izvestiya vuzov. Problemy energetiki, (1-2), 27-36.
8. Shigapov, A. B., Girfanov, A. A., & Shirmanov, M. V. (2012). Calculation of radiative heat fluxes to the walls of the furnace in P5 - approximation of the method of spherical harmonics. Vestnik Kazanskogo gosudarstvennogo tekhnicheskogo universiteta im. A. N. Tupoleva, (1), 18-23.
9. Kuleshov, O. Yu., & Sedelkin, V. M. (2012). Mathematical modeling of local resulting heat transfer in shielded furnaces. Teplovye protsessy v tekhnike, 4(3), 118-124.
Работа поступила Принята к публикации
в редакцию 24.04.2018 г. 28.04.2018 г.
Ссылка для цитирования:
Садыков А. В. К расчету лучистых тепловых потоков в трехмерной прямоугольной области методом дискретных ординат // Бюллетень науки и практики. 2018. Т. 4. №5. С. 1424. Режим доступа: http://www.bulletennauki.com/sadykov (дата обращения 15.05.2018).
Cite as (APA):
Sadykov, A. (2018). To calculation of radiation heat fluxes in three-dimensional rectangular enclosures by method of discrete ordinates. Bulletin of Science and Practice, 4(5), 14-24.
