CC BY
114
©А.В. Садыков, М.А. Бутяков УДК 536.3:535.34
К РЕШЕНИЮ УРАВНЕНИЯ ПЕРЕНОСА ИЗЛУЧЕНИЯ МЕТОДОМ ДИСКРЕТНЫХ ОРДИНАТ
А.В. Садыков, М.А. Бутяков
Нижнекамский химико-технологический институт (филиал) Казанского национального исследовательского технологического университета,
г. Казань, Россия
sadykov_av@mail.ru
Резюме: В работе метод дискретных ординат применяется для решения задачи переноса излучения в трехмерной постановке. Получены результаты с использованием S4-, S6-, S8-приближений и проведено их сравнение с данными других авторов.
Ключевые слова: интенсивность излучения, радиационный теплообмен, контрольный объем.
TO THE SOLUTION OF THE EQUATION OF THE RADIATION TRANSFER TRANSFORMATION BY THE METHOD OF DISCRETE ORDINATES
A.V. Sadykov, M.A. Butyakov
Nizhnekamsk Chemical Technology Institute (branch), Kazan National Research Technological University, Kazan, Russia
sadykov_av@mail.ru
Abstract: In work the method of discrete ordinates is applied to the solution of a problem of transfer of radiation in three-dimensional statement. Results with use of S4-, S6-, S8-approximations are received and comparison with data of other authors is carried out them.
Keywords: intensity of radiation, radiation heat exchange, control volume.
Тепловые потоки в радиантных камерах технологических трубчатых печей формируются в результате совместного действия всех механизмов переноса энергии: лучистого, конвективного и кондуктивного. При температурах, характерных для трубчатых печей нефтехимической промышленности, лучистый теплообмен является основным. Поэтому точность теплового расчета таких установок в основном определяется корректностью модели лучистого переноса энергии. Расчет лучистого теплообмена сводится к решению интегро-дифференциального уравнения переноса излучения (УПИ). В задачах лучистого теплообмена часто используется метод дискретных ординат (МДО). Первоначальный вариант МДО, как развитие двухпотокового приближения Шустера и Шварцшильда для плоской геометрии, был предложен Виком [1] и в дальнейшем развит Чандрасекаром [2]. МДО получил развитие в работах зарубежных и отечественных авторов. Он удачно объединяется с алгоритмами, основанными на применении контрольных
объемов. Метод был использован автором статьи в ряде работ [3 - 5 и др.]. Подробное описание методики расчета для двумерной задачи приведено в статье [5].
Рассмотрим УПИ для области, показанной на рис.1. Для излучающей, поглощающей и рассеивающей серой среды уравнение имеет вид
8/ (Г, О) (г, О) 8/ (г, О)_
8х
8у
ц-
81
= а/ь (Т) - (а + р)/(г, О) +
-Л
4 %
I/(г, О'ЖО', .
(1)
(4%)
Здесь ц, 2,, П - направляющие косинусы; 1(г, О) - интенсивность изучения,
зависящая от положения и направления; а, Р - коэффициенты поглощения и рассеяния
среды соответственно; у(О ', О) - индикатриса рассеяния; /ь (Т) - интенсивность излучения
черного тела при температуре Т .
Уравнение (1) дополняется граничными условиями. Граничное условие к уравнению (1) с учетом диффузного излучения и отражения имеет вид
/(г, О) = г/ь (Т) + Р | \п О'/ (г, О ')сЮ'
(2)
' п а <о
для таких направлений О , для которых п • О > 0. Здесь О ' - направление падающего
излучения; О - направление испускаемого излучения; 8 - степень черноты граничной поверхности; р - отражательная способность поверхности; п - единичный вектор внутренней нормали к границе.
\
У
Рис.1. Система координат
В МДО угловое распределение интенсивности излучения аппроксимируется значениями вдоль ограниченного количества выделенных направлений Б {т=1, 2, ...,
N0} в каждой точке пространства. Для трехмерного поля излучения эти направления задаются набором угловых координат (цт, Ьт, ц„; m= 1, 2, ..., N0}, равных величине проекции единичного вектора направления на оси координат 0х, 0у, 0z соответственно.
В зависимости от количества выделенных направлений различают ^-приближение (N„=8), ^-приближение (N,=24), ^-приближение (N,=48), ^-приближение (N,=80) и другие. Конечно, лучшими приближениями являются те, которые дают точные результаты и не требуют слишком большого объема вычислений. Интегро-дифференциальное уравнение (1) приближенно заменяется системой дифференциальных уравнений относительно интенсивности излучения вдоль каждого из этих направлений:
т=1, 2, 3,...Д0 ,
где - угловые весовые коэффициенты. Индексы т', т обозначают направления
падающего и испускаемого излучения соответственно. Угловые весовые коэффициенты ~№т связывают между собой интенсивности излучения вдоль различных направлений. Эти коэффициенты численно равны площади на поверхности единичной сферы, отсекаемой
соответствующим направлению телесным углом Д^. Коэффициенты и
направляющие косинусы удовлетворяют условиям:
т=1
Индикатриса рассеяния, характеризующая рассеяние лучистой энергии мельчайшими частицами сажи во всех направлениях, записывается в упрощенном виде:
Коэффициент фт,т, учитывающий анизотропию рассеяния при квадратурном представлении интегрального члена, определяется по выражению
Показателем асимметрии является коэффициент g1, который заключен в диапазоне
1 <1. Для изотропного рассеяния gl = 0 (фЫт =1).
Граничное условие (2), характеризующее излучение стенок и отражение лучей, падающих со всех направлений, для различных стенок аппроксимируется следующим образом:
Лт Лт Лт о
Мт + ^т + Чт = Ъ1Ь - (а + Р)/т + ^ £ ^т'Фт'т^т' , (3) дх ду дж 4л т»=!
у(9) = 1 + ^ сов 9.
Фт'т = 1 + ^1[МтМт'+^ т
при х = 0 для мт' < 0 и Мт > 0, при х = Ьх для Мт' > 0 и цт< 0;
О
К = £1ъ {К ) + - Е т К'
Я т'=1
при у = 0 для Е,т> < 0 и £т > 0, при у = Ьу для Е,т> > 0 и £т < 0;
Р
1 т = £1Ъ ) + " Е №т' |Лт ' \{т'
к т=1
при г = 0 для пт' < 0 и пт > 0, при г = Ьг для пт' > 0 и пт < 0.
(6)
В выражениях (4)-(6) вторым членом в правых частях описывается отраженный поток лучистой энергии, при этом суммирование ведется только по направлениям падения лучей.
Область интегрирования покроем конечно-разностной сеткой. Уравнение (3) проинтегрируем по контрольному объему, показанному на рис. 2. В результате получаем следующее конечно-разностное уравнение относительно интенсивностей излучения в узловых точках:
ъ Та
Мтп Ар (1т !т ) + £тВр (1т 1т ) + ЛтСр (1т 1т ) = Ер V р^ + . (7)
р V р^т + °р
Здесь Ар, Вр, Ср, V р - известные функции координат и оптических свойств среды
в точке Р расчетной области Б; Ер - источниковый член в точке Р, зависящий от интенсивности собственного излучения среды; 8р - источниковый член, определяющий
вклад в интенсивность излучения в данном направлении вследствие рассеяния по другим направлениям.
9 Ь
а
с
Рис.2. Контрольный объем 28
Решение системы уравнений (7) совместно с граничными условиями находится итерационным методом. В каждой итерации используется метод покоординатной прогонки. Для получения дополнительных уравнений предположим, что интенсивность излучения в центральной узловой точке контрольного объема связана со значениями в соседних узловых точках следующими выражениями:
^т = ют+(1 - ®) 1Ст = ют+(1 - = ют+(1 - «^т, (8)
где Ю - интерполяционный коэффициент. Формулы для решения алгебраических уравнений получаются следующим образом. Рассмотрим направление интегрирования, для которого |хт > 0, Ьт > 0, пт > 0. Из уравнения (7) с помощью соотношений (8) исключим интенсивности в узловых точках Ь, с, /. В результате приходим к следующей прогоночной формуле:
Тп МтАр 1т + р1 т + ЧтСр1т + ю(( + Бр )
т р + МтАр + ^тБр + Ч тСр
Значения интенсивности в узловых точках Ь, с, / затем вычисляются с помощью соотношений (8). Для других значений угловых координат цт, Ьт, Пт получаются аналогичные выражения.
Итерационный процесс для решения системы уравнений (7) совместно с дискретными аналогами граничных условий (4), (5), (6) реализуется по следующей схеме:
1) задается начальное приближение для интенсивности излучения 1т в узловых
точках расчетной области Б;
2) в граничных узловых точках области производится расчет интенсивностей излучения по выражениям (4), (5), (6). В первой итерации границы считаются черными, а слагаемые, описывающие рассеяние излучения во внутренние ячейки, равными нулю. В последующих итерациях используются заданные значения степени черноты стенок и учитываются члены, описывающие рассеяние излучения на частицах сажи;
3) для каждого ординатного направления (цт, Ь,т, Пт) (т=1, 2, 3, ..., N0) обходятся
все узлы разностной сетки и вычисляются значения интенсивностей 1т по соответствующей прогоночной формуле с учетом знаков цт, Ьт, Пт ;
4) полученные значения интенсивностей 1т принимаются за начальные значения
для следующего шага итерации и происходит переход к пункту 2. Итерационный процесс продолжается до выполнения условия
тах
ре£>
Фр -Фр
фр
<5,
N0
где п - номер предыдущей итерации; 5 - заданная малая величина; ф = £ -
ур
т=1
объемная плотность энергии излучения в точке Р.
После завершения итерационного процесса компоненты вектора плотности результирующего потока излучения в направлениях х, у, z в узловых точках расчетной области определяются суммированием по всем направлениям:
N0 N0 N0 qx = Z ; qy = Z Smwm7m; qz = Z Smwm7i ■
m=1 m=1 m=1
УПИ решается совместно с уравнением энергии, в котором конвекция и теплопроводность предполагаются пренебрежимо малыми:
divq^ = S,
где S - источник тепла. Дивергенция потока излучения находится интегрированием по
всем направлениям Q . Итерационный процесс реализуется по следующему алгоритму:
1) задается начальное поле температуры;
2) решается УПИ при известном поле температуры;
3) решением уравнения энергии находится новое поле температуры.
Пункты 2, 3 повторяются до достижения сходимости по температуре.
Расчеты проведены для идеализированной печи с размерами: Lx =2 м, Ly=2 м, Lz=4 м; данные для среды а = 0,5 м-1 ; S = 5,0 кВт/м3; данные для границ: z=0: Tw=1200 K, е=0,85; z=Lz: Tw=400 K, 8=0,70; другие границы: Tw=900 K, 8=0,70. Здесь Tw - температура соответствующей стенки.
Угловые координаты и весовые множители такие же, как в работе [6]. В расчетах использовалась сетка 7*7*11. Интерполяционный коэффициент ю=0,5. Численные исследования показали, что дальнейшее дробление сетки в случае однородной среды не приводит к уточнению результатов. Вычислительная схема реализована на Фортране.
На рис. 3 показаны распределения температур, полученные с использованием S4-, S6-, S8- приближений, вдоль прямой y=1,0 м на трех плоскостях z = const и их сравнение с P3-приближением по методу сферических гармоник (МСГ) и зональным методом Хоттеля [7]. Вблизи горячей стенки (z=0,4 м) S4-, S6--приближения менее точны, чем P3- приближение МСГ. В центре печи (z=2,0 м) решение, полученное в S6 - приближении, отличается от зонального решения в пределах 1% . Вблизи холодной стенки (z=3,6 м) максимальное отличие решений, полученных в S4-, S6-, S8- приближениях, от зонального решения не больше 1,0%. Следует заметить, что отклонения уменьшаются с ростом оптической толщины. Решения, полученные в S4-, S6-, S8- приближениях отличаются от зонального решения в пределах 0,5% для той же области при коэффициенте поглощения среды а = 1,0 м-1. Результаты расчетов согласуются с численными экспериментами работы [6].
Распределения тепловых потоков на стенках z=0 м и z=4 м вдоль прямой y=1,0 м, полученные с использованием S4-, S6-, S8- приближений, показаны на рис. 4. При а = 0,5 м-1 результирующие потоки излучения, полученные в S4-, S6- приближениях, отличаются от точного решения [7] в пределах 3,0% на горячей стенке (z=0 м) ив пределах 4,0% на холодной стенке (z=4 м). Результаты, полученные в S8 - приближении, не превышают 2,0% на обеих стенках. Результаты расчетов по МДО находятся в лучшем согласии с точным решением, чем по P3- приближению МСГ.
Преимуществом МДО, как уже было отмечено, является сочетаемость с алгоритмами, основанными на применении контрольных объемов. Кроме этого важным преимуществом МДО является то, что для перехода от низшего приближения к более высокому достаточно изменить в расчетной программе значение N0, массивы цт, Ът, nm, wm (m=1, 2, 3, ..., N0). Применяемые для расчета радиационного теплообмена в топках
МСГ [8, 9] и зональный метод [10] более трудоемки. Кроме того, в нерассеивающих средах МСГ дает некорректные результаты [8].
О 0.5 1:0 1.5 х, м 1
Рис. 3. Распределения температур вдоль прямойу=1,0 м на плоскостях г=0,4; г=2,0; г=3,6 м.
Рис. 4. Распределения тепловых потоков на горячей и холодной стенках
Таким образом, исследования показывают, что Б4- приближение МДО дает приемлемые результаты. Для более точных расчетов рекомендуется использовать £6- или £8- приближение. Достигнуто согласование результатов расчетов с численными экспериментами других авторов.
Литература
1. Wick G. Über ebene Diffusionsprobleme // Zeitschrift für Physik. 1943. Bd.121. Pp.702-718.
2. Чандрасекар С. Перенос лучистой энергии. М.: ИЛ, 1953. 431с.
3. Садыков А.В. Методика расчета сопряженного теплообмена в трубчатой печи производства водорода в рамках дифференциального метода // Известия вузов. Проблемы энергетики. 2013. №7-8. С.3-11.
4. Садыков А.В., Вафин Д.Б. Неравномерности обогрева реакционных труб и распределений температуры продуктов сгорания по глубине технологической трубчатой печи // Тепловые процессы в технике. 2014. №8. С.349-355.
5. Садыков А.В. К расчету лучистых тепловых потоков в прямоугольных областях методом дискретных ординат // Известия вузов. Проблемы энергетики. 2016. №3-4. С.13-21.
6. Fiveland W. A. Three-Dimensional Radiative Heat Transfer Solutions by the Discrete-Ordinates Method // Journal of Thermophysics and Heat Transfer. 1988. Vol.2, No. 4. Pp. 309-316.
7. Menguc M., Viskanta R. Radiative Transfer in Three Dimensional Rectangular Enclosures // Journal of Quantium Spectroscopy and Radiative Transfer. Vol. 33, 1985. Pp. 533-549.
8. Шигапов А.Б. Радиационный теплообмен в топках котлов // Известия вузов. Проблемы энергетики. 2014. №1-2. С.27-36.
9. Шигапов А.Б., Гирфанов А.А., Ширманов М.В. Расчет радиационных тепловых потоков к стенкам топки котла в ^-приближении метода сферических гармоник // Вестник Казанского государственного технического университета им. А.Н. Туполева. 2012. № 1. С. 18-23.
10. Кулешов О.Ю., Седелкин В.М. Математическое моделирование локального результирующего теплообмена в экранированных топках // Тепловые процессы в технике. 2012. Т.4, №3. С.118-124.
Авторы публикации
Садыков Айдар Вагизович - канд. техн. наук, доцент, декан факультета управления и автоматизации Нижнекамского химико-технологического института (филиал) «КНИТУ». E-mail: sadykov_av@mail.ru.
Бутяков Марат Анатольевич - аспирант кафедры ЭТЭОП Нижнекамского химико-технологического института (филиал) «КНИТУ».
References
1. Wick G. Über ebene Diffusionsprobleme // Zeitschrift für Physik. 1943. Bd.121. Pp.702-718.
2. Chandrasekhar S. Transfer of radiant energy. Moscow: IL, 1953. 431p.
3. Sadykov A.V. A method for calculating the conjugate heat transfer in a tube furnace for the production of hydrogen in the framework of a differential method, Izvestiya Vuzov. Problems of energy. 2013. № 7-8. C.3-11.
4. Sadykov AV, Vafin D.B. Uneven heating of reaction tubes and temperature distributions of combustion products along the depth of the technological tube furnace. // Thermal processes in engineering. 2014. № 8. P.349-355.
5. Sadykov A.V. To the calculation of radiant heat fluxes in rectangular regions by the method of discrete ordinates // Izvestiya Vuzov. Problems of energy. 2016. № 3-4. C.13-21.
6. Fiveland W. A. Three-Dimensional Radiative Heat Transfer Solutions by the Discrete-Ordinates Method // Journal of Thermophysics and Heat Transfer. 1988 Vol. 4. Pp. 309-316.
7. Menguc M., Viskanta R. Radiative Transfer in Three Dimensional Rectangular Enclosures, Journal of Quantium Spectroscopy and Radiative Transfer. Vol. 33, 1985. Pp. 533-549.
8. Shigapov A.B. Radiation heat exchange in boiler furnaces // Izvestiya Vuzov. Problems of energy. 2014. №1-2. P.27-36.
9. Shigapov AB, Girfanov AA, Shirmanov MV Calculation of the radiative heat fluxes to the walls of the furnace of the boiler in the P5 approximation of the method of spherical harmonics // Bulletin of the Kazan State Technical University. A.N. Tupolev. 2012. № 1. S. 18-23.
10. Kuleshov O.Yu., Sedelkin V.M. Mathematical modeling of local resultant heat exchange in shielded furnaces. Thermal processes in engineering. 2012. T.4, №3. P.118-124.
Authors of the publication
Sadykov Ajdar Vagizovich - the dean of management and automation faculty, candidate of technical sciences, assiciate professor of the department of mathematics
Butyakov Marat Anatol'evich - post graduate student of the department of ETEOP
Поступила в редакцию 1 апреля 2017 г.
