CC BY
46
Вычислительные технологии
Том 19, № 5, 2014
Использование метода конечных объёмов для решения уравнения переноса теплового излучения
в трёхмерных задачах*
К. Ю. Литвинцев1
Метод конечных объёмов для решения уравнения переноса излучения все чаще используется для численного моделирования теплового излучения. Это связано с простотой вывода метода, его гибкостью и отсутствием существенных ограничений на применение. Использование метода конечных объёмов для трёхмерных пространственных задач на произвольных сетках имеет свои особенности, некоторые из них рассматриваются в представленной статье.
Ключевые слова: перенос теплового излучения, математическое моделирование, метод конечных объёмов.
Введение
Для решения задачи переноса теплового излучения к настоящему времени разработано много методов, которые условно можно разделить на ряд групп. Первая группа методов основана на представлении исходного интегро-дифференциального уравнения переноса излучения (УПИ) в форме дифференциальных уравнений второго порядка — Pn-приближения метода сферических гармоник. Наиболее простой из этих методов — Р1-приближение (диффузионный метод), которое может быть получено в том числе прямым интегрированием УПИ по всему телесному углу [1].
В качестве второй группы методов решения уравнения радиационного переноса можно выделить методы Монте-Карло [2]. В общем случае это численный метод решения математических задач при помощи моделирования случайных величин [3]. При использовании данных методов применительно к решению задач переноса лучистой энергии распространение излучения представляется в виде случайной марковской цепи столкновений фотонов или дискретных порций энергии с веществом, приводящих либо к его поглощению, либо к рассеянию.
Третья группа — зональные методы. Основой этих методов является деление излучающей объёмной среды и ограничивающих поверхностей на ряд зон, теплофизические параметры в пределах каждой из которых считаются постоянными [4]. Каждая из зон связана с другими через обобщённые угловые коэффициенты излучения, показывающие, какая доля энергии, излученной в одной зоне, достигает другую зону и поглощается ею.
Четвёртая группа — методы дискретных направлений, в которых проводится не только пространственная, но и угловая дискретизация УПИ. В этих методах всё уг-
1 Институт теплофизики им. С.С. Кутателадзе СО РАН, Новосибирск, Россия
Контактный e-mail: sttupick@yandex.ru
* Работа выполнена при финансовой поддержке Междисциплинарного интеграционного проекта СО РАН на 2012-2014 годы № 49.
ловое пространство разбивается на ряд дискретных телесных углов, внутри которых интенсивность излучения постоянна. К данной группе примыкают потоковые методы, метод дискретных ординат (discrete-ordinates method, DOM) и метод конечных объёмов (finite-volume method, FVM) [5, 6].
Отдельно можно выделить разработанный Локвудом и Шахом [7] метод дискретного переноса, который является гибридным подходом, имеющим многие черты потоковых методов и методов Монте-Карло.
В существующих сегодня наиболее распространённых программных комплексах, таких как ANSYS, Star-CD/CCM+, OpenFAM, FDS, основанных на решении уравнений вычислительной гидродинамики, используются в основном четыре метода решения УПИ: Р1-приближение, FVM, DOM и метод дискретных направлений.
В настоящей работе для решения уравнения переноса теплового излучения рассматривается метод конечных объёмов, основные положения которого впервые были представлены в работе G.D. Raithby, E.H. Chui в 1990 г. [8]. Метод конечных объёмов для решения УПИ был реализован автором представленной статьи в рамках программного комплекса SigmaFlow [9]. Данный метод является оптимальным для моделирования переноса теплового излучения. Так, FVM в отличие от Pn-приближений (особенно Р1) и метода дискретных направлений не имеет ограничений по оптическим свойствам среды (неоднородность, анизотропность, рассеяние), а в отличие от потоковых методов не предъявляет специальных требований к расчётной области. По точности FVM несколько уступает лучевым методам (например методам Монте-Карло), однако существенно выигрывает по требованию к вычислительным ресурсам. При сравнении с дискретно-ординатным методом FVM является немного более ресурсоёмким, но вместе с тем корректнее описывает рассеяние излучения. Ещё одна важная особенность конечно-объёмного метода — возможность использования угловой многоблочности, которая позволяет менять угловую дискретизацию по пространству, что особенно важно для задач с сильной неоднородностью оптических свойств среды.
1. Уравнение теплового переноса излучения
Уравнение радиационного переноса описывает баланс энергии при данной частоте излучения, поступающего вдоль направления s в малый элемент объёма поглощающей, испускающей и рассеивающей среды [10, 11]:
dIvds s) = -ßv (r) Iv (r, s) + Kv (r) Ibv (r) + 1 (r, s') Ф^ (s', s) dQ', (1)
4n
где Iv — интенсивность излучения; I^v — интенсивность излучения абсолютно чёрного тела (АЧТ); ßv — коэффициент затухания (ßv = av + Kv); av — коэффициент рассеяния; Kv — коэффициент поглощения; Q — телесный угол; r — радиус-вектор; s — угловое направление; Фv — индикатриса рассеяния. Первое слагаемое в правой части уравнения (1) отвечает за ослабление интенсивности излучения в направлении s за счёт поглощения и рассеивания радиационной энергии, второе — за усиление интенсивности излучения за счёт излучения среды, третье показывает вклад радиационной энергии с других направлений s вследствие рассеяния излучения.
Далее индекс v — зависимость функций от частоты, будет опускаться, так как уравнения переноса излучения в случаях, когда учитывается спектральная зависимость ра-
диационных свойств среды или используется приближение серого газа, имеют одинаковый вид.
Плотность потока энергии излучения определяется следующим образом:
Е (г) = у I (г, 8) dn. (2)
4п
Источниковый член в уравнении сохранения энергии за счёт радиационного теплообмена для газовой фазы (дивергенция теплового потока излучения) имеет вид
уд = к (Еь (г) - Е (г)), (3)
где ЕЬ — плотность потока энергии излучения АЧТ.
Простейший пример граничных условий, наиболее часто применяемый для прикладных задач, — диффузионное излучение непрозрачной поверхности. В этом случае интенсивность излучения, покидающего диффузионно испускающую и отражающую поверхность, однородно распределяется по всем направлениям:
I (Г, 8) = £ (г) 1Ь (г) + ^ I I (Г, 87) ^ • П| dnf. (4)
б' •п<0
2. Метод конечного объёма
2.1. Вывод дискретного уравнения переноса излучения для ЕУЫ
Дискретный аналог УПИ для метода конечных элементов получается интегрированием уравнения (1) по контрольному объёму и контрольному углу (рис. 1) [5]. В общем случае пространственная сетка может быть неортогональной и неструктурированной [12-14].
Направление распространения луча описывается в наиболее универсальной декартовой системе координат. Таким образом, угловое направление может быть представлено через единичный вектор:
8 = (вт 9 сов ф) ех + (вт 9 вт ф) еу + (сов 9) ег
(5)
а
Рис. 1. Пространственная и угловая дискретизация уравнения радиационного переноса: а — контрольный угол, б — контрольный объём и ориентация контрольного угла
Для удобства запишем уравнение (1) в следующем виде:
¿I (г,б)
¿8
где Б (г, б) — источниковый член
-в (г) I (г, б) + Б (г, б) ,
Б (г, б) = к (г) 1Ь (г, б) + ^ [ I (г, б') Ф (б', б) ¿П'.
4п I
4п
Интегрируя уравнение (6) по контрольному объёму и контрольному углу и заменяя полученный интеграл по контрольному объёму на интеграл по площади граней, получаем
У У 11 (б1 ■ и) ¿A¿n = J J (-в? + Б1) ¿V¿n, (8)
ДП1 ДА ДП1 ДУ
где и и ДА — вектор нормали и площадь грани контрольного объёма соответственно. Левая часть уравнения (8) описывает входящий и выходящий потоки радиационной энергии через грани контрольного объёма, правая часть отражает поглощение и излучение энергии внутри контрольного объёма и контрольного телесного угла. При дискретизации уравнения (8) в методе конечных объёмов величина интенсивности излучения внутри контрольного телесного угла и контрольного объёма принимается постоянной:
£ сАпь i (б1 ■ ипь) ¿П = (-в11 + Б1) Д VДП1, (9)
пь дп
где Б1 = к1ь +--^ I1' Ф1'1 ДП1', Ф1'1 — усреднённая функция рассеивания из контроль-
4п 1'=1
ного угла I' в контрольный угол /; подстрочный индекс пЬ — грань контрольного объёма с узловой точкой Р. Общее число граней зависит как от размерности задачи, так и от типа пространственной сетки. Так, например, для триангулярных пространственных сеток четыре грани соотнесены с узловой точкой Р.
Окончательная дискретизация уравнения (9) для ГУМ для произвольных сеток может быть записана как
£ 4Аь^ь = (-в!р + Б1) Д, (10)
пЬ
где В1пЬ = J (б1 ■ ипь) ¿П. Для связи между интенсивностями на гранях контрольной ДП1
ячейки с интенсивностью внутри контрольного объёма используется схема "по потоку", при которой значение интенсивности в контрольном объёме сносится на грани, расположенные по ходу распространения луча в направлении б1 . Дискретная форма ГУМ для схемы аппроксимации "по потоку" для произвольных сеток записывается следующим образом:
ь
а1Р Рр = £ аПь4 + Ь1, (11)
пь
где
alP = £ max (AnbDlnb, 0) + вр AVPАП1,
nb
alnb = max (-A,nbDnb, 0) , bl = SPAVpAQl.
Кроме схемы простой аппроксимации "по потоку" можно использовать экспоненциальную схему [5, 13]. Суть данной аппроксимации заключается в связи интенсивности внутри контрольного объёма с интенсивностью на грани вдоль распространения излучения согласно закону Бугера — Ламберта — Бера.
Помимо различных схем аппроксимации для FVM существуют различные способы дискретизации углового пространства. Простейший способ — однородная угловая дискретизация (см. рис. 1, б). Более сложная угловая дискретизация предполагает неравномерное разделение углового пространства по азимутальному и полярному углам для достижения более равномерного распределения величины дискретного телесного угла. Так, в одной из модификаций FVM — FTn в зависимости от значения полярного угла предполагается увеличивать или уменьшать кратно двум количество угловых разбиений вдоль азимутального угла [15].
2.2. Граничные условия
2.2.1. Диффузионное излучение непрозрачной поверхности
Диффузионное излучение непрозрачной поверхности предполагает, что излучение с поверхности носит изотропный характер как для отражённого, так и для испускаемого поверхностью излучения и описывается уравнением (5), дискретный аналог которого для FVM имеет вид [5]
iw = е ■ Ib + Р £ IW ■ DW, (12)
fc,D4 <0
где индекс l означает исходящее от поверхности, индекс k — падающее излучение.
2.2.2. Периодические граничные условия
При периодических граничных условиях связываются две границы, на которых задаётся условие равенства распределения входящей, для одной границы, и исходящей — для другой, интенсивности излучения. Сложность реализации периодических граничных условий для FVM заключается в возникновении несоответствия контрольных углов на периодических границах при их повороте относительно друг друга (рис. 2) [16]. Для решения этой проблемы вводится локальная система координат, получающаяся из глобальной путём применения матрицы трансформации для перехода от одной периодической границы к другой. Для локальной системы строится угловая дискретизация с таким же разбиением, как и для глобальной. Далее для одной из систем все телесные углы разбиваются на малые углы, с помощью которых определяется область пересечения телесных углов из разных координатных систем. На основании этого определяется связь между телесными углами на периодических границах / :
Dl > 0, Il = £ /'I1', £//' = 1, i' i'
Периодические границы
/ I \\У Контрольный угол,
" АО.1* возникающий при повороте г, периодической границы
Существующий контрольный угол
Рис. 2. Несоответствие контрольных углов при повороте периодической границы
в' > о, I1' = £ аI1, = 1.
(14)
Здесь индексы I и I' означают направления распространения излучения для интенсив-ностей на соответствующих периодических границах.
2.3. Метод "бегущего" счёта
Для нахождения поля интенсивностей для уравнений (11) используется метод "бегущего" счёта [5, 13], суть которого заключается в последовательном расчёте интенсивности излучения в контрольных объёмах в выбранном направлении. Для любой сетки можно построить такую последовательность обхода контрольных объёмов для произвольного направления распространения излучения, что излучение, поступающее в данный контрольный объём, будет всегда известно. При использовании этого метода на декартовых сетках при условии абсолютно чёрных стенок и отсутствия рассеяния в среде поле интенсивности находится за один итерационный шаг.
3. Особенности использования метода конечных объёмов для пространственных задач переноса излучения
В настоящей статье не рассматриваются проблемы, связанные с лучевым эффектом, возникающим при недостаточной дискретизации углового пространства, или с выбором разбиения углового пространства. Исследуются только два аспекта работы с ГУМ: ошибки, возникающие вследствие использования произвольных пространственных сеток при малой угловой дискретизации, и использование угловой многоблочности.
3.1. Перекрытие контрольных углов с гранями контрольных объёмов
В результате использования произвольных сеток возникает проблема перекрытия контрольных углов с гранями контрольных объёмов, когда границы последних не совмещены с прямоугольной системой координат (рис. 3) [5, 16]. При пересечении телесного угла гранью контрольного объёма лучистый поток разбивается на две части: входящий
и исходящий. При использовании системы уравнений (12), (13) значения интенсивности для входящего П+ и исходящего П- лучистых потоков принимаются одинаковыми: либо 1р, либо 1пь. В реальности значения интенсивности исходящего и входящего потоков разные, что в случае сильно анизотропных сред может приводить к существенным ошибкам.
3.1.1. Локальное расщепление телесного угла
Для уменьшения ошибки, вносимой пересечением граней контрольных объёмов с контрольными телесными углами, телесный угол П1 предлагается разбивать плоскостью грани на два угла П+ и П-, соответствующие входящему и выходящему лучистым потокам. В этом случае уравнение (11) модифицируется следующим образом:
апЬ = — Апб^П+,
= / (8' • Ппь) ^
д«г±
а1Р = £ АпьО П- + вр АУр АП1,
пЬ
Ь1 = ^ АУр АП1. (15)
Здесь надстрочные индексы + и — соответствуют входящему в контрольный объём и исходящему из него радиационному потоку.
Рис. 3. Пример перекрытия контрольного угла с гранью контрольного объёма
3.1.2. Влияние процедуры локального расщепления телесного угла на точность решения
Рассмотрим влияние перекрытия контрольных углов с гранями контрольных объёмов и процедуры локального расщепления телесного угла на результаты расчёта поля излучения ГУМ на примере тестовой задачи.
Тестовая задача представляет собой замкнутую прямоугольную область с различной постоянной температурой стенок: 1400 К — нижняя, 400 К — верхняя, 273 К — боковые стенки. Поле температуры рассчитывается решением уравнения сохранения энергии. Степень черноты холодной стенки полагается равной единице, для остальных стенок — 0.7. Принимаются два значения коэффициента поглощения: к = 1.0 и к = 0.1 м-1, при этом критерий Бугера составляет соответственно 1.6 и 0.16. Телесное пространство для ГУМ разбивается на 32 контрольных угла. Для оценки влияния эффекта углового перекрытия и процедуры локального расщепления телесного угла рассматриваются два способа расположения объекта относительно глобальной декартовой системы координат: в первом случае грани контрольных объёмов совпадают с ней, эффект углового перекрытия отсутствует, во втором — объект повернут на 30° относительно оси Оу и возникает эффект углового перекрытия (рис. 4). В рамках задачи проводились три варианта расчёта: первый — базовый, когда угловое перекрытие отсутствует, второй — когда присутствует угловое перекрытие, но не используется процедура локального расщепления телесного угла, в третьем варианте присутствуют и угловое перекрытие, и процедура локального расщепления телесного угла.
В целом при повороте объекта и соответственно появлении углового перекрытия наблюдаются смещение максимума, нарушение симметрии и изменение величины поля плотности потока энергии относительно базового варианта (рис. 5). Однако отклонения решения от базового варианта для различных значений коэффициента поглощения сильно отличаются. Так, при к = 1.0 м-1, когда среда является близкой к оптически толстой, симметрия нарушается слабо, но в среднем величина плотности потока излуче-
Рис. 4. Геометрия и сетка объекта (сетка 9 х 9 х 19)
90000 —,
80000 —
70000 —
60000
И | I | I
0.4
0.8 1.2
X, м
п I 1 I 1.6 2
Рис. 5. Распределение плотности потока энергии излучения вдоль линии (см. рис. 4): а — к = 1,
б — к = 0.1 м-1; • — без углового перекрытия, А — с коррекцией углового перекрытия, +--
без коррекцией углового перекрытия
а
ния в отсутствие коррекции на 27 % превышает базовый вариант, а при использовании коррекции соответственно не превышает 3% (рис. 5, а). При к = 0.1 м-1, когда среда является оптически прозрачной и возникает сильная анизотропия поля излучения, наблюдается значительное нарушение симметрии, но при этом уровень отклонения величины плотности потока излучения уменьшается. Так, без использования коррекции отклонение составляет 4.8%, а с использованием — 2.4% (рис. 5, б). Сильное нарушение симметрии при к = 0.1 м-1 указывает на недостаточную угловую дискретизацию, с увеличением числа контрольных углов этот эффект уменьшается.
Таким образом, применение процедуры локального расщепления телесного угла позволяет не только повысить точность расчёта, но и уменьшить его время.
При появлении эффекта углового перекрытия количество внутренних итераций, необходимых для построения поля излучения в случае фиксированных значений полей остальных физических величин, при решении УПИ по сравнению с базовым вариантом возрастает. Увеличение числа внутренних итераций с использованием и без использования процедуры коррекции происходит по-разному (см. таблицу). При этом в случае отсутствия локального расщепления телесных углов существует зависимость количества внутренних итераций от оптических свойств среды: для к = 1.0 м-1 по сравнению с к = 0.1 м-1 относительное время возрастает почти в 1.5 раза. При использовании процедуры коррекции такая зависимость не наблюдается: относительное время в обоих случаях одно и то же.
Сравнение времён расчёта относительно базового варианта
Значение к, м 1 Базовый вариант Без коррекции С коррекцией
1.0 1.0 2.05 1.3
0.1 1.0 1.39 1.3
На рис. 6 представлены графики сходимости плотности потока энергии излучения по внешним итерациям (16) при изменении полей физических величин для к = 1.0 и к = 0.1 м-1:
е> =
Е (Е - Е-1)
1=1,ы
Е Е
(16)
где г — номер контрольного объёма (И — число контрольных объёмов), 3 — внешний номер итерации.
Графики сходимости базового варианта и варианта с процедурой коррекции для различных оптических сред близки (см. рис. 6). Если коррекцию не применять, то при к = 1.0 м-1 наблюдается уменьшение числа внешних итераций до достижения сходимости. Последнее связано с тем, что в данном случае число внутренних итераций сильно возрастает. Так, если ограничить время окончания счёта не сходимостью (см. таблицу), а количеством итераций, то время расчёта относительно базового варианта увеличивается в 2.6 раза. Полученные данные указывают на то, что использование процедуры коррекции как по точности решения, так и по времени расчёта является наиболее эффективным для оптически более толстой среды. Процедура локального расщепления телесного угла, как правило, наиболее важна в случае применения малой угловой дискретизации и существенной оптической анизотропии среды. Таким образом, использование данной процедуры позволяет в некоторых случаях отказаться от избыточной дискретизации пространства и тем самым существенно уменьшить вычислительные затраты на решение УПИ.
Рис. 6. Графики сходимости плотности потока энергии излучения для к = 1 м-1 (а) и к = 0.1 м-1 (б): 1 — без углового перекрытия, 2 — с угловым перекрытием без процедуры коррекции, 3 — с угловым перекрытием с процедурой коррекции
3.2. Угловая многоблочность
3.2.1. Описание связи блоков на границах
Метод конечных объёмов даёт возможность использования угловой многоблочности, позволяющей задавать разную угловую дискретизацию по рассматриваемому пространству [17].
Стыковка блоков на границе контрольных объёмов с разной угловой дискретизацией происходит из условия сохранения потока энергии излучения:
Я
' ДА ^ 2п
I (8 • п) dШA
блок1
'ДА ./2п
I (8 • п) dШA
:17)
блок2
Из уравнения (17) выводится связь интенсивностей на общей грани контрольных объёмов с различной угловой дискретизацией:
ДA,
пЬ
V Т1 В1
пЬ пЬ
1=1,Дг , >0
' по
ДА,
пЬ
блок1
к
Е
_Тк вк - пЬ пЬ
к=1,ЭК <0
по
:18)
блок2
Таким образом, связь интенсивности в направлениях I и к для различных блоков можно описать следующим образом:
V ^ а 1к в
/ у Т■nЬJ■nЬВ
пЬ пЬ пЬ
г=1,д',- >о
Г_ Iк Вк 1
пЬ пЬ\ блок2'
:19)
блок1
где ¡г — доля энергии излучения, уходящей из телесного угла ДП1 и попадающей в ДПк через грань пЬ:
¡г
1к
в1к
ВпЬ
В1 .
ВпЬ
(20)
Применение угловой многоблочности позволяет сокращать время расчёта в случае предварительно известных оптических свойств среды. Так, если в задаче присутствуют расчётные области, требующие хорошей угловой дискретизации, то без использования угловой многоблочности детальную угловую дискретизацию необходимо распространять на всю расчётную область, что приводит к существенному увеличению времени расчёта.
Ниже рассмотрена тестовая задача, в которой оптические свойства среды неоднородны и для каждой из областей можно использовать собственную угловую дискретизацию.
3.2.2. Демонстрация использования угловой многоблочности
Расчётная область в тестовой задаче имеет кубическую форму (рис. 7, а), среда в ней принимается неизлучающей, т.е. к = 0.0м-1. В области х < 0.5м среда считается изотропно рассеивающей: коэффициент рассеивания а полагается равным 2 м-1, в остальной части куба а = 0 м-1. На нижней стенке задаётся поток излучения д, верхняя и боковые стенки считаются абсолютно чёрными и неизлучающими (см. рис. 7, а) [14].
Представлены результаты расчётов трёх вариантов: первый — 32 телесных угла, второй — 288 телесных углов, третий — 32 телесных угла в области а = 2 м-1 и 288 —
0.24 —I
0.2-
0.16-
О}
0.12-
0.08-
0.04«
• 288 направлений
---32
- 32+228
"Г
0.2
"Г
0.4
0.6
Т
0.8
а
Рис. 7. Пример использования угловой многоблочности: а — постановка задачи, б — распределение безразмерного падающего потока излучения вдоль линии г = 1 м, у = 0.5 м
в области а = 0 м-1. При использовании малого количества угловых разбиений появляется "нефизичный" минимум падающего потока в области без рассеивания, где в эталонном решении с 288 телесными углами наблюдался максимум, средняя величина ошибки достигает 11.5%, а максимальная — 30% (см. рис. 7, б), но при этом время расчёта уменьшается в 36 раз. Непропорциональный уменьшению числа контрольных углов рост скорости расчёта связан с наличием источникового члена, описывающего энергию излучения, поступающую за счёт рассеяния с других направлений излучения. При использовании процедуры многоблочности удаётся снизить отклонение от эталонного решения в среднем до 4.6 %. Наибольшее расхождение наблюдается на границе сред с различными значениями коэффициента рассеяния. Временные затраты на расчёт по сравнению с эталонным решением уменьшились примерно в 3 раза, несмотря на увеличение на 17% числа внешних итераций, необходимых для достижения заданной сходимости. Таким образом, процедура угловой многоблочности позволяет существенно сократить расчётное время решения УПИ без существенной потери точности.
Заключение
В работе показано, что достоинства и недостатки метода конечных объёмов определяются аппроксимацией углового пространства: с одной стороны — простота вывода, гибкость при разбиении углового пространства, с другой — сложности работы с произвольными пространственными сетками, лучевой эффект и ошибки моделирования рассеивания [5, 19]. При этом достоинства ГУМ можно использовать для его модификаций как для снижения вычислительных затрат, так и для существенного уменьшения влияние ошибок угловой дискретизации. На представленном примере показано, что в случае возникновения пространственного перекрытия контрольных углов и объёмов, управляя
угловой дискретизацией, можно уменьшить эту ошибку без значительного увеличения числа контрольных углов. Другой рассмотренный в статье инструмент управления угловой дискретизацией, существующий в FVM, — угловая многоблочность, которая, локально задавая угловую дискретизацию, позволяет учитывать пространственное распределение оптических свойств среды.
Таким образом, для решения УПИ метод конечных объёмов является не только эффективным, но и гибким, позволяющим адаптировать его применительно к условиям рассматриваемой задачи [5, 6, 18].
Список литературы
[1] Четверушкин Б.Н. Математическое моделирование задач динамики излучающего газа. М.: Наука, 1988. 304 с.
[2] Snegirev A.Y. Statistical modeling of thermal radiation transfer in buoyant turbulent diffusion flames // Intern. J. of Heat and Mass Transfer. 2004. Vol. 136. P. 51-71.
[3] Соволь И.М. Метод Монте-Карло. М.: Наука, 1988. 64 с.
[4] Блох А.Г., Журавлёв Ю.А., Рыжков Л.Н. Теплообмен излучением: Справочник. М.: Энергоатомиздат, 1991. 432 с.
[5] Chai J.C., Patankar S.V. Finite-volume method for radiation heat transfer // Advances in Numer. Heat Transfer. Taylor and Francis Group, 2000. Vol. 2. P. 109-138.
[6] Litvintsev K.Yu., Dekterev A.A. Comparison of the finite-volume and discrete-ordinate methods and diffusion approximation for the radiative heat transfer equation // Heat Transfer Res. 2008. Vol. 68. P. 653-655.
[7] Lockwood F.C., Shah N.G. New radiation solution method for incorporation in general combustion prediction procedures // Proc. of Eighteenth Symp. (Intern.) on Combustion. Pittsburg, 1981. P. 1405-1414.
[8] Raithby G.D., Chui E.H. A Finite-volume method for predicting a radiant heat transfer enclosures with participating media // J. of Heat Transfer. 1990. Vol. 11. P. 415-423.
[9] Дектерёв А.А., Гаврилов А.А., Минаков А.В. Современные возможности CFD кода SigmaFlow для решения теплофизических задач // Современная наука: Исследования, идеи, результаты, технологии. 2010. Вып. 2(4). C. 117-122.
[10] Modest M.F. Radiative Heat Transfer. Second Edition. London, New York: Acad. Press. Elsevier Sci., 2003. 842 p.
[11] Howell J.R., Siegel R., Menguc M.P. Radiative Heat Transfer. Thermal Radiation Heat Transfer 5th edn. New York: CRC Press. Taylor and Francis Group, 2010. 957 p.
[12] Murthy J.Y., Mathur S.R. Finite volume method for radiation heat transfer using unstructured meshes //J. of Thermophys. and Heat Transfer. 1998. Vol. 12, No. 3. P. 313-321.
[13] Chai J.C., Rath P. Discrete-ordinates and finite-volume methods for radiation heat transfer // Intern. Workshop on Discrete-Ordinates and Finite-Volume Methods for Radiation Heat Transfer. Guwahati, 2006. URL: http://hdl.handle.net/2080/341
[14] Литвинцев К.Ю. Особенности использования конечно-объёмного, дискретно-ординат-ного методов и диффузионного приближения для уравнения радиационного теплопере-носа // Вестник Сибирского гос. аэрокосм. ун-та. 2008. № 4(21). C. 44-47.
[15] Guedri K., Naceur B.M., Rachid M., Rachid S. Formulation and testing of the FTn finite volume method for radiation in 3-D complex inhomogeneous participating media // J. of Quant. Spectroscopy and Radiative Transfer. 2006. Vol. 98. P. 425-445.
[16] Moder J.P., Kumar G.N., Chai J.C. An unstructured-grid radiative heat transfer module for the national combustion code // AIAA 38th Aerospace Sciences Meeting & Exhibit. Reno. 2000. AIAA-2000-0453.
[17] Chai J.C., Moder J.P. Radiation heat transfer calculation using an anglur-multiblock procedure // Numer. Heat Transfer. 2006. Vol. 38. P. 1-13.
[18] Ramamoorthy B., Cheng G.C., Koomullil R.P., Rahmani R.K. Finite volume method for non-equilibrium radiative heat transfer // Intern. J. of Heat and Mass Transfer. 2013. Vol. 65. P. 670-681.
[19] Tan H., Zhang H., Zhen B. Estimation of ray effect and false scattering in approximate solution method for thermal radiative transfer equation // Numer. Heat Transfer. 2004. Part A, vol. 46. P. 807-829.
Поступила в редакцию 15 апреля 2014 г., с доработки — 14 августа 2014 г.
Application of the finite-volume method for 3D radiation heat transfer problems
Litvintsev Kirill Yu.1
Purpose. Implementation of the finite-volume method for 3D radiation heat transfer problems has some peculiarities that are addressed in the present research. Control-angle overlap and angular-multiblock procedure is of particular interest.
Methodology. Control-angle overlap is encountered when edges of the control volumes cross control angles. It leads to errors and increased calculation time when the angle of discretization is coarse. Errors appear as a result of the neglected difference between incoming and outgoing radiant energy for an overlapping control angle. Two principal approaches can help to solve this problem: to refine the angle of discretization or to subdivide the "bad" control angles. This article introduces the second approach when an overlapping control angle is divided into two sub-control angles by the edge of the control volume. Angular-multiblock procedure allows changing angular discretization in accordance with optical properties of the environment thus calculation time is reduced. Originality/value. Solution of the test problem, which addresses the overlap control angles with the faces of control volumes, has shown that the division into two parts of the control angle reduces the computational error, as well as significantly reduce the computation time. To test the procedure we reduve the initial problem to the corner multiblock problem in which the scattering coefficient varied greatly in space. For this task, three alternatives were considered angular rate: 288, 32, 288 + 32. The first option — the reference solution. For the second variant, we obtained a nonphysical solution. In the third case, a solution was obtained close to the reference, and the calculation time was reduced in three times.
Findings. The finite-volume method is flexible. It allows extension of its application scope
to computations of the radiation heat transfer processes.
Keywords: radiative heat transfer, numerical simulation, finite-volume method.
Received 15 April 2014 Received in revised form 14 August 2014
1 Institute of Thermophysics SB RAS, 630090, Novosibirsk, Russia Corresponding author: e-mail: sttupick@yandex.ru
Aknowlegements: Interdisciplinary Integration Project of SB RAS During 2012-2014 years, No. 49.
