CC BY
35
УДК 536.3
К РЕШЕНИЮ УРАВНЕНИЯ ПЕРЕНОСА ИЗЛУЧЕНИЯ МЕТОДОМ ДИСКРЕТНЫХ ОРДИНАТ В ОСЕСИММЕТРИЧНОЙ ЦИЛИНДРИЧЕСКОЙ ОБЛАСТИ
А.В.Садыков
Нижнекамский химико-технологический институт (филиал) Казанского национального исследовательского технологического университета, г. Нижнекамск, Россия, sadykov_av@mail.ru
Резюме: В работе метод дискретных ординат применяется для решения задачи переноса излучения в конечной осесимметричной цилиндрической области. Получены решения для условий, моделирующих поглощающую среду. Приведены результаты, полученные с использованием S2-, S.f-приближений, и проведено их сравнение с экспериментальными данными других авторов.
Ключевые слова: интенсивность излучения, радиационный теплообмен, метод дискретных ординат, контрольный объем.
TO THE SOLUTION OF THE EQUATION OF TRANSFER OF RADIATION BY METHOD OF DISCRETE ORDINATES IN THE AXISYMMETRIC DOMAIN
A.V.Sadykov
Nizhnekamsk Institute of Chemical Technology (Branch) of " Kazan National Research
Technological University " Nizhnekamsk , Russia, sadykov_av@mail.ru
Abstract: In this work the method of discrete ordinates is applied for the solution of the problem of transfering of radiation in a finite axisymmetric domain. Decisions for the conditions to model an absorbing medium are obtained. The results received with S2-, S4- approximations are given and their comparison with the experimental data of other scientists is carried out.
Keywords: intensity of radiation, radiation heat exchange, discrete ordinates method, control volume.
В нефтехимической промышленности используются цилиндрические трубчатые печи для нагрева сырья реактора. При условиях, характерных для трубчатых печей, радиационный теплообмен является основным. При заданном поле температуры, при известных радиационных свойствах среды и ограничивающих поверхностей расчет лучистого теплообмена в стационарном случае сводится к решению интегро-дифференциального уравнения переноса излучения (УПИ) при соответствующих граничных условиях. Современные методы расчета переноса излучения описаны в работе [1]. В задачах лучистого теплообмена часто используется метод дискретных ординат (МДО). Метод привлекателен тем, что он удачно объединяется с алгоритмами, основанными на применении контрольных объемов. В УПИ для криволинейной геометрии, в отличие от прямоугольной, присутствуют производные по угловым переменным и возникает проблема приближенной оценки этих производных. Имеются различные способы аппроксимации этих производных. Наиболее удачным является способ [2]. В указанной работе уравнение переноса нейтронов записывается в дивергентной форме, и разностные
уравнения выводятся методом баланса частиц в ячейке. Дивергентная форма приводит к конечно-разностным схемам более высокого качества [2, 3].
УПИ для излучающей, поглощающей и рассеивающей серой среды для цилиндрической геометрии в дивергентной форме имеет вид
И • djrl(M, П)) -1. 8(ц1 (M, П)) 8I(M, П) = ^ TT) -
r dr r 8y 8z
(a + ß)I(M, П)+A f I(M, Q)y(Q', n)dQ'. (1)
4n d) _
Здесь ц, , n - направляющие косинусы; I(M, П) - интенсивность излучения в точке M в
направлении a, ß — коэффициенты поглощения и рассеяния среды соответственно;
у(П', П) — индикатриса рассеяния; (T) — интенсивность излучения черного тела при
температуре T. Система координат показана на рис. 1. Для осесимметричной цилиндрической геометрии интенсивность излучения в любой точке M излучающего объема зависит от двух пространственных координат r, z и от двух угловых переменных 0 и I(M, П) = I(r, z, 0, ¥).
V г
Рис. 1. Система координат для цилиндрической печи. Обозначения: 1 - трубчатый экран; 2 - боковая
стена; 3 - горелка
Уравнение (1) дополняется граничными условиями. Граничное условие к уравнению (1) на жесткой стенке с учетом диффузного излучения и отражения имеет вид
I(М, П) = г1ь (Т) + Р | \п -Й'I(М, ПУ П' (2)
% ПП'<0
для таких направлений О , для которых п ■ О > 0 . Здесь О' - направление падающего
излучения; О - направление испускаемого излучения; 8 - степень черноты граничной поверхности; р - отражательная способность поверхности; п - единичный вектор внутренней нормали к границе. Такое граничное условие ставится при г = Я, г = 0, г = Ь. На оси симметрии (г = 0) ставится условие симметричности решения.
Индикатриса рассеяния, характеризующая рассеяние лучистой энергии мельчайшими частицами сажи во всех направлениях, записывается в следующем виде:
у(О', О) = 1 + ^СОБ ф.
В МДО угловое распределение интенсивности излучения аппроксимируется
значениями вдоль ограниченного количества выделенных направлений От {т = 1, 2, ...,
N0} в каждой точке пространства. Эти направления задаются набором угловых координат
(цт, Ът, пт; т = 1, 2, ..., N0}. Каждое направление От изображается точкой на
поверхности единичной сферы. Количество направлений в ^„-приближении равно N0=„(„+2), где „ - порядок приближения. Интеграл в уравнении (1) приближенно заменяется квадратурной формулой, а уравнение (1) - системой связанных дифференциальных уравнений относительно интенсивности излучения вдоль каждого из
направлений От:
^т • д(Пт) 1 д(Цт!т) +, д1 т = . , , Р) . + Р Д м ф .
- • -Д---• -I- +4т Т--Шь - (а + р)1т +— Е Мт'фт'тУт' , (3)
г дг г ду дг 4к т=1
т=1,2,3,.,N0 ,
где Мт - угловые весовые коэффициенты; цт, пт - направляющие косинусы для
дискретного направления От: цт = sin0cosY, пт = sin0sinY, Ьт = cos0. Индексы т', т
обозначают направления падающего и испускаемого излучения соответственно. Угловые весовые коэффициенты wm связывают между собой интенсивности излучения вдоль различных направлений. Эти коэффициенты численно равны площади на поверхности
единичной сферы, отсекаемой соответствующим направлению От телесным углом ДО.
Коэффициент фт 'т, учитывающий анизотропию рассеяния при квадратурном
представлении интегрального члена, определяется выражением
фт 'т :
Показателем асимметрии является коэффициент g1, который заключен в диапазоне - 1 ^ <1. Для изотропного рассеяния g1 = 0 (фт 'т =1).
Граничное условие (2) для различных стенок аппроксимируется следующим образом:
р ^0
4 = 81Ь (Тм ) + ~ Е Мт' ^т' 1 т' (4)
Л т '=1
при г = Я для цт' > 0 и цт< 0;
(т =81Ь (Тм ) + Р
при г = 0 для Е,т ' < 0 и Е,т > 0;
1т = 81 ь (Тм ) + РЕ Мт т'К' (5)
Л т '=1
Р Щ
1т =г1Ь Т ) + "Е ' £т' 4 '
К т=1 (6)
при г = Ь для > 0 и < 0.
В выражениях (4)-(6) вторым членом в правых частях описывается отраженный поток лучистой энергии, при этом суммирование ведется только по направлениям падения лучей. Граничное условие при г =0:
1т = 1т , ^т=Цт'. (7)
В уравнении (3) присутствует производная по угловой переменной Используя способ, предложенный в работе [2] для аппроксимации такой производной, второй член в левой части уравнения (3) представляется в виде
1. дСЛю^ю) = 1 , т , , т , )
----ит+1/21ш+12 ит-1/27ю_12) ' wm ,
г г
где аю±^2 - неизвестные коэффициенты.
Область интегрирования покроем конечно-разностной сеткой. Обе части уравнения (3) умножим на 2п г d г d г и проинтегрируем по контрольному объему. В результате получим следующее конечно-разностное уравнение относительно интенсивностей излучения в узловых точках:
^т (А+1Л+1 _ А1/1 ) _ (А +1 _ А )(ат+1/2Im+V2 _а'т_1/2^-1/2 ) / ^т + ^
+4т (В у+! Iу+! - Ву/у ) = ЫУ1Ь _ (« + Р)К/т + ^^ .
Здесь Аг-, В у - площади поверхностей граней ячейки в радиальном и осевом направлениях; V - объем ячейки; 8р - источниковый член. Для изотропного рассеяния
ву N0 8 = ву Е ^ / ,.
р л ' 1 т т
4п ю'=1
Рекуррентные соотношения для нахождения неизвестных коэффициентов
выводятся из условия сохранения баланса. Цель введения этих коэффициентов состоит в том, чтобы баланс сохранялся и для конечных интервалов. В уравнении (8) индексы г, т относятся, соответственно, переменным г, г,
Для получения дополнительных уравнений предположим, что интенсивность излучения в центральной узловой точке контрольного объема связана со значениями в соседних узловых точках следующими выражениями:
= ю/+1 + (1 _ = ю/т+1/2 + О _ Ю)/ю_12 = ю/у +1 + (1 _ ю)/у , (9) где ю - интерполяционный коэффициент. При ю=0,5 получаем «ромбовидную» схему, а при ю=1 - «шаговую» схему или схему «против потока».
Формулы для решения алгебраических уравнений получаются следующим образом. Допустим, что цт>0, Е,т > 0. Из уравнения (8) с помощью выражений (9) выразим интенсивность излучения в центральной точке ячейки:
_ УтА11 + Рт/т_У2 + ^тВ// + (УЬ + 8р )
т = + Рт +^тВ + (<х + Р)У '
где А = 4+1 ; В = В/ + Ву+1 ; Рт = _(4+1 _ А )(«т+1/2 + «т-1/2 ) / ^т . Для
других комбинаций ц,т и получаются аналогичные выражения.
Решение системы уравнений (8) совместно с граничными условиями находится итерационным методом. В каждой итерации используется метод покоординатной прогонки.
Для осесимметричной цилиндрической геометрии количество ординатных направлений уменьшается в 2 раза. Сначала уравнение (3) решается в направлении пт=0, и полученное решение принимается за начальное приближение. Вычисления начинаются в правом верхнем углу, используя граничное условие (4) . Для каждого ординатного направления обходятся все узлы сетки и вычисляются значения интенсивностей по соответствующей прогоночной формуле с учетом знаков направляющих косинусов. В зависимости от их знаков счет начинается от разных границ. Изменение знака направляющих косинусов означает перемену соответствующего направления интегрирования. Итерационный процесс продолжается до выполнения критерия сходимости.
«Ромбовидная» схема иногда приводит к отрицательной интенсивности, что невозможно по физическим соображениям. Этого можно избежать, если использовать комбинацию «ромбовидной» схемы и разности «против потока». В процессе вычислений, когда значения интенсивности становятся меньше некоторой заданной величины, включается схема «против потока». При использовании такой схемы появление отрицательной интенсивности исключено.
В расчетах использованы приведенные в [4] угловые координаты и весовые множители, которые заимствованы из работы [5].
Для проверки адекватности математической модели и программного кода проведены тестовые расчеты. Расчеты проведены для цилиндрической печи [6] со следующими размерами: длина 1=5 м, радиус Л=0,45 м. Данные для среды а = 0,3 м-1 ; данные для границ: 7^=425 К, £„,=0,8. В работе [6] имеются экспериментальные данные плотностей тепловых потоков к стенке для незакрученного пламени природного газа. Распределение температуры внутри печи известно [6]. При известном поле температуры и других известных данных решается УПИ.
На рис. 2 показаны распределения плотностей тепловых потоков дя к стенке. Из рис. 2 видно, что результаты ^-приближения находятся в хорошем соответствии с экспериментальными данными [6]. ^-приближение дает чуть завышенные значения (около 4%) плотностей тепловых потоков на участке от 1 до 2,5 м, но при дальнейшем увеличении г дает практически такие же результаты как ^-приближение. ^-приближение, в котором уравнение решается в 24 направлениях в осесимметричной цилиндрической геометрии, не обеспечило улучшения результатов по сравнению с ^-приближением.
Як, -т-■-т-т-
кВт
0 -'-1-'-'-
0 1 2 з г, м
Рис. 2. Распределения плотностей тепловых потоков к стенке
МДО был использован автором статьи в ряде работ [7 - 9 и др.] при решении задачи сложного теплообмена в технологических трубчатых печах цилиндрического и коробчатого типов. В работе [10] МДО использован при моделировании работы акустических горелок в трубчатых печах. Результаты анализа влияния методов решения уравнения переноса излучения на расчетные тепловые характеристики цилиндрической печи типа ВА-301 приведены в статье [9]. В указанной работе процессы, протекающие в топке, моделируются двумерными уравнениями энергии, переноса излучения, движения, k-e модели турбулентности и модели горения природного газа. Для решения уравнения переноса излучения использованы метод сферических гармоник (^-приближение) и метод дискретных ординат (S2-, ^-приближения). На основе сравнения расчетных значений температуры топочных газов с экспериментальными данными сделан вывод о возможности использования этих приближений МДО при расчете тепловых характеристик цилиндрической печи.
Таким образом, S2- и S.i-приближения МДО моделируют радиационный теплообмен в осесимметричной цилиндрической геометрии с приемлемой точностью. Для такой геометрии нет необходимости использования более высоких приближений МДО.
Литература
1. Dombrovsky L.A. The use of Transport Approximation and diffusion - based Models in Radiative Transfer Calculations // Computational Thermal Sciences. 2012. Vol.4. No.4. Pp. 297-315.
2. Carlson B.G., Lathrop K.D. Transport Theory - The Method of Discrete Ordinates: In Computing Methods in Reactor Physics. Edited by H. Greenspan, C.N. Kelber, D. Okrent. Gordon and Breach Science Publishers. New York. 1968.
3. Басс Л.П. Конечно-разностные методы решения уравнения переноса в задачах со сложной геометрией. ИПМ АН СССР. Препринт №14. М.: 1974. 75с.
4. Jamaluddin A.S., Smith P.J. Predicting Radiative Transfer in Axisymmetric Cylindrical Enclosures Using the Discrete Ordinates Method // Combustion Science and Technology. 1988. Vol.62. No.4-6. Pp.173186.
5. Truelove J.S. Evaluation of a Multi-Flux Model for Radiative Heat Transfer in Cylindrical Furnaces. AERE R-9100,AERE Harwell, U.K. 1978.
6. Wu H.L., Flicker N. An Investigation of the Behavior of Swirling Jet Flames in a Narrow Cylindrical Furnace. 2nd' Members Conference, International Flame Research Foundation, Ijmuiden, The Niderlands. 1971.
7. Садыков А.В., Садыкова Д.А., Вафин Д.Б. Влияние параметров горелок и их расположения на аэродинамику топочных газов и тепловые характеристики цилиндрических трубчатых печей // Известия вузов. Проблемы энергетики. 2012. №5-6. С.17-24.
8. Садыков А.В., Вафин Д.Б. Неравномерности обогрева реакционных труб и распределений температуры продуктов сгорания по глубине технологической трубчатой печи // Тепловые процессы в технике. 2014. №8. С.349-355.
9. Садыков А.В. Расчет тепловых характеристик цилиндрической печи с использованием разных методов решения уравнения переноса излучения // Известия вузов. Проблемы энергетики. 2016. № 9-10. С. 43-47.
10. Вафин Д.Б., Бутяков М.А. Трехмерное моделирование работы акустических горелок в трубчатых печах // Известия вузов. Проблемы энергетики. 2016. № 9-10. С. 48-55.
Автор публикации
Садыков Айдар Вагизович - канд. техн. наук, доцент кафедры математики, декан факультета управления и автоматизации Нижнекамского химико-технологического института (филиал) Казанского национального исследовательского технологического университета (КНИТУ). E-mail: sadykov_av@mail.ru.
References
1. Dombrovsky L.A. The use of Transport Approximation and diffusion - based Models in Radiative Transfer Calculations // Computational Thermal Sciences. 2012. Vol.4. No.4. Pp. 297-315.
2. Carlson B.G., Lathrop K.D. Transport Theory - The Method of Discrete Ordinates: In Computing Methods in Reactor Physics. Edited by H. Greenspan, C.N. Kelber, D. Okrent. Gordon and Breach Science Publishers. New York. 1968.
3. Bass L.P. Of course - difference methods of solution of the equation of transfer in tasks with the composite geometry. The Keldysh Institute of Applied Mathematics, Ac. Of Sc., the USSR. The preprints No. 14, Moscow. 1974.P.75.
4. Jamaluddin A.S., Smith P.J. Predicting Radiative Transfer in Axisymmetric Cylindrical Enclosures Using the Discrete Ordinates Method // Combustion Science and Technology. 1988. Vol.62. No.4-6. Pp.173186.
5. Truelove J.S. Evaluation of a Multi-Flux Model for Radiative Heat Transfer in Cylindrical Furnaces. AERE R-9100,AERE Harwell, U.K. 1978.
6. Wu H.L., Fricker N. An Investigation of the Behavior of Swirling Jet Flames in a Narrow Cylindrical Furnace. 2nd. Members Conference, International Flame Research Foundation, Ijmuiden, The Niderlands. 1971.
7. Sadykov A.V., Sadykova D.A., Vafin D.B. Effect of the parameters of burner and their location on aerodynamics of combustion gases and thermal characteristics of cylindrical tubular furnaces // News of HIGH SCHOOLS. Power problems. 2012. № 5-6. P.17-24.
8. Sadykov A.V., Vafin D.B. Unevenness of heating of reactionary pipes and distributions of temperature of products of combustion on depth of the technological tubular furnace // Thermal Processes in engineering. 2014. №8. P.349-355.
9. Sadykov A.V. Calculation of thermal characteristics of the cylindrical furnace with use of different methods of the solution of the equation of transfer of radiation // News of HIGH SCHOOLS. Power problems. 2016. № 9-10. P.43-47.
10. Vafin D.B., Butyakov M.A. Three-dimensional modeling work acoustic burners tubular furnaces // News of HIGH SCHOOLS. Power problems. 2016. № 9-10. P.48-55.
Authors of the publication
Aydar V. Sadykov - cand. sci. (techn.), assiciate professor, Department of mathematics, Dean of management and automation faculty, Nizhnekamsk Institute of Chemical Technology (Branch) of " Kazan National Research Technological University ".
Поступила в редакцию 04 июня 2017 г.
