Использование схемы “бегущего счета” для решения СЛАУ, полученных на основе метода конечных объемов для уравнения теплового переноса излучения Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 22, № 4, 2017

Использование схемы "бегущего счета" для решения СЛАУ, полученных на основе метода конечных объемов для уравнения теплового переноса излучения

К. Ю. ЛитвинцЕв1,2'*, С. А. Филимонов2 1 Институт теплофизики СО РАН, Новосибирск, Россия 2Сибирский федеральный университет, Красноярск, Россия *Контактный e-mail: [email protected]

Метод конечных объемов получил широкое распространение при численном моделировании процессов переноса теплового излучения. Это связано с простотой вывода метода, его гибкостью и отсутствием существенных ограничений на его применение. Работа посвящена возможности использования схемы "бегущего счета" для решения системы линейных алгебраических уравнений, полученной на основе метода конечных объемов для уравнения переноса излучения на неструктурированных сетках, и сравнению его с другими методами, используемыми в вычислительной гидродинамике.

Ключевые слова: перенос излучения, математическое моделирование, метод конечных объемов, схема "бегущего счета".

Введение

Несмотря на то что метод конечных объемов (МКО) для решения уравнения переноса излучения (УПИ) начал интенсивно развиваться только в 90-х годах прошлого века [1, 2], к настоящему времени он получил достаточно широкое распространение. Особенно можно отметить его использование в свободно распространяемых программных продуктах ОрепЕОАМ и ЕВБ, предназначенных для численного моделирования процессов, связанных с гидрогазодинамикой и сложным тепломассопереносом. Метод конечных объемов реализован в авторском программном комплексе SigmaFlow, разработанном в Красноярском филиале ИТ СО РАН и предназначенном для численного исследования и анализа процессов гидрогазодинамики и теплообмена на основе методов вычислительной гидродинамики. Расчеты с использованием МКО, результаты которых представлены в работе, проведены на базе программы SigmaFlow [3, 4].

Одной из основных задач, с которой приходится сталкиваться при использовании МКО, является поиск оптимальных алгоритмов для минимизации вычислительных ресурсов, направленных на решение задачи переноса излучения, поскольку время, затраченное на ее решение, может многократно превышать время, необходимое для численного моделирования остальных процессов. Подходы, направленные на уменьшение времени расчета УПИ методом конечных объемов, основаны на оптимизации углового разбиения телесного пространства, использовании пространственных прямоугольных

© ИВТ СО РАН, 2017

гексагональных сеток, разработке и совершенствовании различных схем аппроксимации, применении методов параллельных вычислений и подборе методов решения систем линейных алгебраических уравнений (СЛАУ) [5-8]. Данная работа посвящена возможности использования схемы "бегущего счета" при решении СЛАУ для неструктурированных многогранных пространственных сеток и сравнению ее с методами решения СЛАУ, применяемыми в вычислительной гидродинамике.

1. Краткое описание метода конечных объемов для решения уравнения переноса излучения

Стационарное уравнение переноса излучения описывает баланс энергии, поступающей вдоль направления в в малый элемент объема поглощающей, испускающей и рассеивающей среды:

^8) = -Р (г) I (г, в) + к (г) 1Ь (г) + ^ [ I (г, в') Ф (в', в) ¿П', (1)

аз J

где I — интенсивность излучения; 1ъ — интенсивность излучения абсолютно черного тела (АЧТ); [3 = к + а — коэффициент затухания; к — коэффициент поглощения; а — коэффициент рассеяния; П — телесный угол; г — радиус-вектор; в — угловое направление; Ф — функция рассеяния.

Для простоты уравнение (1) записано для случая, когда оптические свойства среды не зависят от частоты излучения. Плотность потока энергии излучения определяется следующим образом:

Е (г) = JI (г, в) ¿П. (2)

В качестве граничных условий может использоваться диффузионное излучение непрозрачной поверхности. В этом случае интенсивность излучения, покидающего диффузионно испускающую и отражающую поверхность, однородно распределяется по всем направлениям:

I (г, в) = е (г) 1Ь (г) + ^ J I (г, в') |в' ■ п| ¿П'. (3)

При использовании метода конечных объемов дискретный аналог УПИ получается интегрированием уравнения (1) по контрольному объему и контрольному углу [1, 2]. В итоге дискретный вид уравнения (1) для произвольных пространственных сеток при 1-м направлении распространения излучения может быть записан следующим образом:

£ 11пЬАпЬ^пЬ = (-^ + Б1) ДУДП, (4)

пЬ

где

Мп „

^ /Ф^ДП'; = (в1 ■ ппЬ) dП;

1/_1 ^

^ = к1ъ + ф"ДП; °пъ = I (в1

1'=1 ДПI

X

Рис. 1. Разбиение телесного пространства при Пл-модификации МКО

ДО} — контрольный телесный угол; подстрочный индекс пЬ — грань контрольного объема с узловой точкой Р; Апь — площадь грани пЬ; ипь — вектор нормали грани пЬ; Ф11 — усредненная функция рассеивания из контрольного угла I' в контрольный угол I; Nn — число дискретных телесных углов; ДУ — контрольный объем.

При использовании противопоточной схемы аппроксимации значение интенсивности в контрольном объеме сносится на грани, расположенные по ходу распространения луча в текущем направлении. Дискретная форма МКО для противопоточной схемы аппроксимации в случае использования произвольных сеток записывается следующим образом:

а1Р 11Р = £ а1пЬ 11пЬ + Ь1, (5)

пЬ

где

а1пЬ = тах (-АпЪВ1пЬ, 0) , а1Р = ^ тах (АпЪВ1пЬ, 0) + рРДУРДП1, Ь1 = Б1Р ДУРДП1.

пЬ

Для разбиения углового пространства используется ЕТп-модификация МКО (рис. 1), когда число контрольных углов неравномерно распределено по полярному углу, т. е. дискретизация азимутального угла зависит от значения полярного [9]. В этом случае можно получить более равномерную дискретизацию углового пространства по сравнению с равномерным разбиением полярного и азимутального углов. При таком подходе общее число телесных углов можно описать следующим выражением:

Щ/2

п1-1, если N0 — четное число,

^ = < ¿¡/2-1 (6) 2^0 Е п%-1 + ^п(Ме+1)/2, если N0 — нечетное число,

г=1

где N0 — число разбиений по полярному углу; N0 — начальное число разбиений по азимутальному углу (при г = 0); п — множитель (целое число); = г Д9.

2. Метод решения системы СЛАУ для переноса излучения на основе схемы "бегущего счета"

Для решения СЛАУ можно применять как стандартные методы, используемые при решении уравнений гидродинамики и тепломассообмена (ЫСОЯТЛБ, Б1Ьи, СЯ и т.д.), так и схему "бегущего счета". В первом случае преимущества заключаются в том, что используются одни и те же методы для решения как уравнений гидродинамики, так и для уравнения переноса излучения, соответственно сохраняется общая структура представления данных, кроме того, эти методы работают для произвольных сеток. Основное преимущество схемы "бегущего счета" состоит в том, что это прямой метод решения СЛАУ. Его суть заключается в последовательном расчете интенсивности излучения в контрольных объемах в выбранном направлении (рис. 2). При его использовании строится такая последовательность обхода контрольных объемов в заданном направлении распространения излучения, при которой излучение, поступающее в контрольный объем, будет всегда известно. Построение обхода контрольных объемов возможно только для сеток с выпуклыми многогранными ячейками. Другой недостаток схемы "бегущего счета" — это ограниченный набор схем аппроксимации, при котором его можно использовать, фактически это противопоточная схема и некоторые ее модификации (например, экспоненциальная схема) [2].

Наиболее наглядно продемонстрировать преимущество схемы "бегущего счета" можно при решении задачи переноса излучения, когда в качестве граничных условий принимаются стенки со степенью черноты е = 1 и отсутствует рассеяние в среде, в этом случае поле интенсивности излучения находится за один итерационный шаг.

Алгоритм определения порядка обхода контрольных объемов вдоль заданного направления 1-го распространения излучения для схемы "бегущего счета" состоит из следующих действий:

1) маркировка граней контрольного объема, через которые излучение поступает в ячейку и уходит из нее;

2) определение граней, на которых значения интенсивности входящего изучения известны;

3) определение стартовых контрольных объемов, у которых известны все значения интенсивности входящего излучения, и включение их в список обхода;

4) изменение статуса интенсивности излучения на гранях контрольного объема, включенного в список обхода, с "неизвестного" на "известный";

5) определение следующего контрольного объема, у которого известны все значения интенсивности входящего излучения, и переход к п. 4 до тех пор, пока не будут задействованы все контрольные объемы.

У

Рис. 2. Демонстрация обхода контрольных объемов по схеме "бегущего счета"

3. Постановка задачи

С целью определения эффективности схемы "бегущего счета" для УПИ сравнение проводилось с итерационными методами: методом неполной ЬИ-факторизации (Б1Ьи) и стабилизированным методом бисопряженных градиентов (ЫСОЯТЛБ) [10, 11]. Для приближения полученных результатов к реальным задачам рассматривалось решение совместно уравнений переноса излучения и теплопроводности. Расчетная область представляла собой неправильную Т-образную геометрию, в которой две стенки горячие со степенью черноты е, равной 1.0, а остальные — 0.7 (рис. 3). Принимается, что среда не рассеивает излучение. Поглощающие свойства среды рассматривались от полностью

6.0 м

а

1 3.0 м :

£ = 1.0 ъ £ = 0.7

Т = 603 К Т = 303 К

а = 5 Вт/(К • м) а = 10 Вт/(К- м)

1.0 м

2.0 м

Рис. 3. Геометрия и параметры тестовой задачи

Рис. 4. Используемые типы сеток: а — окто; б — гекса; в — тетра

прозрачной до оптически толстой (г = 100). Температура Т изменялось от 300 до 600 К, коэффициент теплопроводности а варьировался от 5 до 10 Вт/(К • м).

Для анализа использовались три вида различных сеток (рис. 4): гексасетка (шестигранные ячейки), тетрасетка (четырехгранные ячейки) и сетка на основе окторазбиения (состоит в основном из шестигранных ячеек, а также пятигранных ячеек на границах и в областях изменения шага дискретизации). Для каждого типа сеток рассматривались два уровня дискретизации: базовый и детальный. Для октосетки детализация была увеличена в 4.5 раза с 1900 ячеек, для гекса — в 8.4 раза с 2400 ячеек, а для тетра — в 3.7 раза с 1500 ячеек. Расчет останавливался при достижении заданной сходимости по полю температуры.

Необходимо отметить, что количество требуемых итераций для достижения заданной сходимости при использовании различных методов было близким в большинстве расчетов. Сходимость поля излучения определялась по изменению поля плотности потока энергии излучения (см. уравнение (2)).

4. Результаты расчетов

Чтобы оценить эффективность решения только задачи переноса излучения, сначала рассматривался вариант, для которого было заморожено поле температур, а изменялись только два параметра: оптическая толщина и степень черноты стенок (степень черноты принималась для всех стенок одинаковой). Для данной задачи использовалась гексагональная сетка (рис. 4). В результате оказалось, что метод Б1ССЯТЛБ и другие методы "крыловского" типа (CR.ES, СО) не могут находить решения СЛАУ с нулевой правой частью, т. е. когда среда оптически прозрачная (коэффициент поглощения равен нулю). В данной задаче это происходило, когда степень черноты принималась равной 0.7 и ниже. Поэтому дальнейшее сравнение схемы "бегущего счета" проводилось в первую очередь с методом В1Ьи.

На рис. 5 приведены результаты расчетов в виде графика зависимости отношения времен счета с использованием метода В1Ьи и схемы "бегущего счета" (¿*) от оптической толщины среды при различных значениях степени черноты стенок. Как видно, максимальная эффективность "бегущего счета" достигается при степени черноты, равной 1.0, что является очевидным результатом, так как в этом случае поле излучения рассчитывается за один проход (одну итерацию). При этом максимальное превосходство над методом ВГЬи до 11 раз достигается в условиях оптически тонкой среды. Когда поглощательные свойства среды усиливаются, ¿* снижается примерно в шесть раз. При задании степени черноты, отличной от 1.0, ¿* в среднем становится меньше и не является уже монотонно убывающей. Это связано с необходимостью итерационного счета для получения поля излучения с использованием схемы "бегущего счета" из-за перераспределения энергии излучения со стенок вследствие появления отражения. Примечательно, что уменьшение степени черноты с 0.7 до 0.5 фактически не привело к изменению зависимости ¿* от оптической толщины среды. Кроме того, на графиках видно, что при увеличении оптической толщины величина ¿* для е = 1.0 и е = 0.7 (0.5) принимает близкие значения.

Далее рассматривалась задача, в которой совместно решались уравнения переноса излучения и теплопроводности (см. рис. 3) с применением различных сеток (см. рис. 4). Как видно из рис. 6, несмотря на использование различных сеток и дискретизации, рассчитанные поля плотности потока излучения для представленных вариантов близки.

Рис. 6. Поле плотности потока излучения при г =10 для трех различных типов сетки: тет-ра (а), гекса (б), окто (в), Вт/м2. Шаг изолиний 1000 Вт/м2

Результаты сравнения времени расчета с использованием метода Б1Ьи и схемы "бегущего счета" при совместном решении уравнений переноса излучения и теплопроводности представлены на рис. 7. Первое, что необходимо отметить, величина £* в среднем

меньше, чем в предыдущем примере (см. рис. 5), когда рассчитывалось только поле излучения. Фактически время расчета во второй задаче определяет сходимость решения уравнения теплопроводности, скорость которой падает с ростом оптической толщины и соответственно увеличением радиационного источникового слагаемого в правой части.

Изменение детализации сетки также влияет на эффективность "бегущего счета" по сравнению с Б1Ьи. Так, если для оптически тонкой среды £* больше для базовых сеток, то для оптически толстой среды наоборот, при этом количество итераций до достижения заданной сходимости практически не зависит ни от типа сетки, ни от ее детализации. Результаты по влиянию угловой дискретизации не приведены ввиду их незначительности, так, например, при уменьшении угловой дискретизации с 640 контрольных углов до 240 это приводит в основном к небольшому увеличению £* (менее 5 %). Если рассматривать в целом результаты, приведенные на рис. 7, то £* с ростом т преимущественно падает, хотя и не монотонным образом, вне зависимости от типа сетки. Это можно объяснить тем, что вычислительные затраты на одну итерацию с методом Б1Ьи падают с ростом оптической толщины, в то время как при использовании "бегущего счета" они относительно слабо колеблются около постоянной величины (рис. 8). Здесь (£)* = • (Ът) (и), где (Ът) — среднее время на одну итерацию для расчетов с заданной

т

оптической толщиной г, Ит — число расчетных вариантов с разными значениями г.

Таким образом, уменьшение вычислительных затрат с использованием Б1Ьи при росте т приводит к уменьшению £*, а "произвольные" колебания затрат времени на одну итерацию для "бегущего счета" обусловливают немонотонность этого уменьшения. Данные "произвольные" колебания появляются вследствие взаимодействия двух разнонаправленных тенденций, возникающих при увеличении оптической толщины среды: пер-

Рис. 7. Сравнение времени расчета с использованием метода Б1Ьи и схемы "бегущего счета" для решения СЛАУ применительно к УПИ. Количество контрольных углов 640

Рис. 8. Относительные средние вычислительные затраты на одну итерацию

вый связан с уменьшением количества внутренних итераций на расчет поля излучения, а второй — с ухудшением сходимости решения уравнения теплопроводности.

Из графика, представленного на рис. 8, видно, что если для схемы "бегущего счета" указанные разнонаправленные тенденции примерно уравновешивают друг друга, то в случае использования Б1Ьи на общее время расчета сильнее влияет уменьшение количества внутренних итераций для поля излучения, чем ухудшение сходимости решения уравнения теплопроводности. Особым случаем является вариант с прозрачной средой, когда во всех случаях максимум £* приходится вблизи т = 0 (см. рис. 5 и 7), т. е. максимальная эффективность схемы "бегущего счета" по сравнению с методом Б1Ьи достигается для оптически тонкой, но непрозрачной среды.

Помимо представленных в статье расчетов проводились вычисления поля излучения на более сложной геометрии с использованием октосетки с малой угловой дискретизацией (24 телесных угла) и сопряженным теплообменом для оптически прозрачной среды, которые также показали, что использование схемы "бегущего счета" позволяет снизить время на решение задачи по сравнению с методом Б1Ьи примерно в три раза.

Сравнение результатов расчетов с использованием методов Б1Ьи и ЫСОЯТЛБ (для оптически непрозрачных сред) показало, что в большинстве случаев расчеты с Б1СО-ЯТЛБ-методом требуют больше времени для достижения решения, для некоторых вариантов более чем в 1.5 раза.

Заключение

В целом схема "бегущего счета" как прямой метод для решения СЛАУ, полученных на основе метода конечных объемов для уравнения переноса излучения, показал свою эффективность по сравнению с более универсальными методами, использующимися

в вычислительной гидродинамике, позволяя в 2.5-5 раза уменьшать время, требуемое на расчет поля излучения на произвольных выпуклых многогранных сетках. Однако необходимо отметить, что эффективность схемы "бегущего счета" будет падать в случае малой угловой дискретизации при использовании процедур уменьшения ошибки углового перекрытия (пересечения гранью контрольного объема контрольного угла) [4]. В этом случае необходимо брать значения интенсивности излучения с "неизвестных" контрольных объемов, и схема становится частично итерационной.

Благодарности. Работа выполнена при финансовой поддержке РФФИ, Правительства Красноярского края, Красноярского краевого фонда поддержки научной и научно-технической деятельности в рамках научного проекта № 16-48-242085.

Список литературы / References

[1] 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. 112. P. 415-423.

[2] 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. USA: American Institute of Aeronautics and Astronautics, 2000. 9 p. AIAA 2000-0453.

[3] Дектерев А.А., Гаврилов А.А., Минаков А.В. Современные возможности CFD кода SigmaFlow для решения теплофизических задач // Современная наука: исследования, идеи, результаты, технологии. 2010. Вып. 2(4). C. 117-122.

Dekterev, А.А., Gavrilov, А.А., Minakov, А.У. New features of SIGMAFLOW code for thermophysics problems solving // Modern Science: Researches, Ideas, Results, Technologies. 2010. Iss. 2(4). P. 117-122. (In Russ.)

[4] Литвинцев К.Ю. Использование метода конечных объемов для решения уравнения переноса теплового излучения в трехмерных задачах // Вычисл. технологии. 2014. Т. 19, № 5. C. 37-50.

Litvintsev, K.Yu. Application of the finite-volume method for 3D radiation heat transfer problems // Comput. Technologies. 2014. Vol. 19, No. 5. P. 37-50. (In Russ.)

[5] 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. 15 p.

[6] Wray, A. Improved finite-volume method for radiative hydrodynamics // Seventh Intern. Conf. on Comput. Fluid Dynamics (ICCFD7). Big Island, Hawaii, July 9-13, 2012. 20 p.

[7] Lygidakis, G.N., Nikolos, I.K. Improving the accuracy of a finite-volume method for computing radiative heat transfer in three-dimensional unstructured meshes // 3rd SouthEast Europ. Conf. on Comput. Mechanicsan ECCOMAS and IACM Special Interest Conf. Kos Island, Greece, 12-14 June, 2013. P. 599-620.

[8] Trovalet, L., Jeandel, G., Coelho, P.J., Asllanaj, F. Modified finite-volume method based on a cell vertex scheme for the solution of radiative transfer problems in complex 3D geometries // J. of Quant. Spectroscopy & Radiative Transfer. 2011. Vol. 112. P. 2661-2675.

[9] 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.

[10] Barrett, R., Berry, M.W., Chan, T.F. et al. Templates for the solution of linear systems: Building blocks for iterative methods. Philadelphia: SIAM, 1994. 107 p.

Available at: http://www.netlib.org/templates/templates.pdf

[11] Van der Vorst, H.A. Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for solution of non-symmetric linear systems // SIAM J. Sci. Stat. Comp. 1992. No. 2. P. 631-644.

Поступила в 'редакцию 9 января 2017 г., с доработки — 1 июня 2017 г.

The use of the "marching order" scheme for solving SLAE obtained on the basis of the finite volume method for the radiating heat transfer equation

LlTVINTSEV, KlRILL Yu.1'2'*, FlLIMONOV, SERGEY A.2

institute of Thermophysics SB RAS, Novosibirsk, 630090, Russia 2 Siberian Federal University, Krasnoyarsk, 660041, Russia * Corresponding author: Litvintsev, Kirill Yu., e-mail: [email protected]

Purpose. An opportunity of using the "marching order" scheme to solve system of linear algebraic equations (SLAE) applied to the finite-volume method for the radiation transfer equation on unstructured grids and a comparison with iterative methods used in CFD are the subject of this research.

Methodology. Both the standard iterative methods utilized in CFD and the "marching order" scheme can be used to solve the SLAE. The iterative methods both can be used to solve the hydrodynamic equations and the radiation transfer equation. They work for arbitrary meshes. The "marching order" scheme is a direct method to solve SLAE, which is its main advantage. The essence of "marching order" is the sequential direct calculation of the radiation intensity in the control volumes in the chosen direction. The "marching order" scheme has been compared with iterative methods to determine its efficiency compared to with the method of incomplete LU-factorization (DILU) and the biconjugate gradient stabilized method (BiCGSTAB).

Originality/value. The "marching order" scheme has turned to be faster 4 to 11 times than the DILU method in the case of solving the radiation transfer task alone. It has been faster from 2 to 6 times compared to iterative methods in case of the combined equations of radiation transfer and heat conduction.

Findings. The "marching order" scheme has shown its efficiency in comparison with more universal iterative methods for solving the SLAE. It allows reducing the calculation time of the radiation field on arbitrary convex polyhedral meshes in 2.5-5 times.

Keywords: radiative heat transfer, numerical simulation, finite-volume method, marching order.

Acknowledgements. This research was supported by RFBR, Government of Krasnoyarsk Territory, Krasnoyarsk Region Science and Technology Support Fund to the research project No. 16-48-242085.

Received 9 January 2017 Received in revised form 1 June 2017

© ICT SB RAS, 2017