Использование графического конвейера для расчета угловых коэффициентов излучения

И.А. Лебедев,

соиск., LebedevMSZ@gmail.com,

ИПУ РАН, г. Москва

В статье описывается новый подход к расчету угловых коэффициентов излучения, необходимых при моделировании тепловых процессов зональными методами. Принципиальной особенностью нового алгоритма является использование графического конвейера для большей части вычислений. Суть метода заключается замене источник излучения «наблюдателем» и оценке угловых размеров объектов системы, попадающих в поле его зрения. Тем самым определяется, какая часть излучения, распространяемая из источника, попадает на эти объекты, то есть определяются угловые коэффициенты излучения.

The article describes a new approach to the calculation of angular radiation coefficients required for modeling of thermal processes zonal methods. The principal feature of the new algorithm is to use a graphics pipeline for most of the calculations. The method is the replacement of the radiation source, "an observer" and evaluation of the angular dimensions of objects falling into his field of vision. This determines which part of the radiation propagating from the source falls on the objects, that is determined by the angular coefficients of radiation.

Введение

Практика показала, что наиболее приемлемым и перспективным с точки зрения расчётно-теоретического анализа работы теплотехнических установок являются зональные методы. Было разработано несколько схем расчета с целью эффективного применение зонального метода к различным теплотехническим задачам. Эти схемы характеризуются различием математической записи выражений, аппроксимирующих интегральные уравнения излучения, последовательностью и способами их решения[1-5]. При определении лучистой составляющей зональных уравнений теплового баланса в нашей стране распространена расчетная схема, использующая введенные Ю.А. Суриновым разрешающие угловые коэффициенты излучения [6].

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

Особенностью системы уравнений с разрешающими угловыми коэффициентами излучения является простота формы записи и универсальность уравнений. Одним из главных недостатков метода является сложность расчета угловых коэффициентов, особенно в случае сложной геометрической формы исследуемого оборудования. Используемый в последнее время метод Монте-Карло при всех своих преимуществах обладает слишком большой вычислительной ёмкостью. Вместо существующих способов расчета предлагается использовать принципиально новый графический способ расчета. Каждый из элементарных элементов системы разбивается на некоторое число локальных источников излучения. Интенсивность излучения каждого отдельного локального источника уменьшается пропорционально их количеству относительно интенсивности излучения всего элемента. В построенной заранее трехмерной модели исследуемого теплового объекта в позицию каждого локального источника излучения ставится камера (наблюдатель). Дальнейшим шагом является отрисовка (рендеринг) средствами DirectX полученного изображения в сферической системе координат (в этом случае обзор будет охватывать всю полусферу излучения). Всё, что попадает в поле зрения камеры, принимается как приёмник излучения источника, а отношение площади, занимаемой приёмником, ко всей площади изображения, принимается равным угловому коэффициенту излучения. Использование стандартной графической библиотеки DirectX 11 фирмы Microsoft для построения изображения позволит переложить большую часть расчетов на графический процессор.

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

Угловой коэффициент излучения – это величина, характеризующая поток излучения, который попадает с одной площадки изучаемой системы на другую. Если i-тая излучающая поверхность (источник излучения) системы излучает на j-тую излучающую поверхность системы (приёмник излучения), то определяет, какая часть от общего излучения i-той поверхности попадает на j-тую. В литературе часто приводится формула (рисунок 1), которая, однако плохо работает при переходе к алгебраическим уравнениям, в случае если площади поверхностей сопоставимы с квадратом расстояния между ними:

где, - углы между нормалями излучающих поверхностей и прямой, соединяющей их центры, - расстояние между излучающими поверхностями.

- площадь поверхности-приёмника.

Рис.1. – К определению углового коэффициента

Если принять излучение полностью диффузионным (независящим от направления), то угловой коэффициент отражает, какая часть излучения источника, распространяемого по полусфере, попадает на проекцию приёмника на эту сферу. Используя термины трехмерной компьютерной графики, получается, что если поставить наблюдателя в точку источника излучения, и направить взгляд по направлению нормали в данной точки и отобразить полученное изображение в сферической системе координат (зенитные и азимутные углы откладываются по осям x и y, r выполняет роль глубины изображения – координаты z), то отношение полученной площади приёмника к полной площади изображения, есть угловой коэффициент излучения. Из этого следует, что для расчета углового коэффициента излучения необходимо:

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

2) Направлением камеры выбрать нормаль к излучающей поверхности в данной точке.

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

4) Определить на полученном изображении отношение площадей видимых поверхностей к площади изображения (точнее круговой центральной его части, так как координаты точек будут иметь разброс по первым двум координатам от -π/2 до π/2).

Для реализации алгоритма предлагается использовать стандартные средства для работы с графическими процессорами, а именно DirectX 11 фирмы Microsoft.

Рис. 2. Графический конвейер

2. Реализация на основе графического конвейера

Графический процессор это отдельное устройство в составе компьютера, которое выполняет обработку графической информации для вывода на экран (рендеринг). Благодаря специализированной конвейерной архитектуре он справляется с этой задачей намного лучше центрального процессора. Видеоадаптеры (видеокарты) имеют в своем составе до нескольких тысяч ядер и собственную оперативную память, их архитектура оптимизирована для массивных параллельных вычислений, что и определяет их большую эффективность в рендеринге. Графический конвейер – это аппаратно-программный комплекс визуализации трёхмерной графики, последовательность стадий (функций), выполняющихся в фиксированном порядке и параллельно. Каждая стадия принимает информацию из предыдущей стадии и отправляет в следующую. Стандартный графический конвейер обрабатывает вершины, геометрические примитивы, а также пиксели конвейерным способом. В отличие от центрального процессора, графический процессор имеет большое количество ядер, что позволяет ему обрабатывать высокий объем данных параллельно, не ставя их в очередь. В современном графическом конвейере поддерживаемом в DirectX 11 пять программируемых стадий (рисунок 2): вершинный шейдер(Vertex Shader), шейдер поверхности(Hull Shader), шейдер домена(Domain Shader), геометрический шейдер(Geometry Shader) и пиксельный шейдер(Pixel Shader).

Каждая из этих стадий необходима для реализации алгоритма. Прежде чем запускать графический конвейер в видеокарты необходимо загрузить текстуры (с точки зрения программирования это обычный массив) и буфер вершин. В буфере трехмерная модель, которую необходимо отрисовать, представлена в виде последовательности точек (вершин), объединённых в треугольники. Получив данные видеокарта запускает графический конвейер, в котором каждая из стадий выполняется параллельно на всех имеющихся ядрах (шейдерных процессорах).

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

Так как графический конвейер способен работать только с вершинами в декартовой системе координат, то он будет некорректно проводить растеризацию, то есть заполнения пространства внутри треугольника. Прямые, соединяющие вершины в декартовой системе координат, не являются таковыми в сферической, вырождаясь в дуги, что вызывает погрешности, вплоть до полного искажения результатов. Решением данной проблемы в рамках графического конвейера может быть увеличение числа треугольников, то есть уменьшение их размеров. Но задавать заранее большее количество вершин накладно с точки зрения ресурсов и отнимает у алгоритма гибкость (необходимо загружать разные модели для разной точности расчетов). Более правильным является использование двух следующих стадий конвейера, появившихся в 11 версии DirectX: шейдер поверхности и шейдер домена. Их совокупная работа позволяет организовать так называемую тесселяцию – разделение поступающих на вход треугольников на более мелкие треугольники «на лету», т.е. генерируемые вершины не присутствуют в исходной модели, а появляются на этой стадии. Это позволяет резко повысить детализацию трехмерной модели, в том числе и в сферической системе координат, и следственно точность расчета угловых коэффициентов. Чем больше уровень тесселяции тем точнее грани первичного треугольника соответствуют их реальному виду в сферической системе координат. Код этих шейдеров стандартен, и в статье не приводится. Данный способ позволяет организовать изменение уровня тесселяции в зависимости от расстояния до объекта и его углового размера, так как нет смысла тратить ресурсы графического процессора на деление далеких от камеры треугольников, они и так отобразятся близко к истине, в то время как близкие объекты будут искажены, и требуют более высокого уровня тесселяции. Так же теселяция позволяет на «на лету» округлять поверхности, что удобно при прорисовке сфер, цилиндров и других сферических объектов, так как они при разбиении на излучающие поверхности заменяются правильными многоугольниками и многогранниками.

Полученные вершины поступают на следующую стадию – геометрический шейдер. В отличие от вершинного он обрабатывает весь примитив целиком, т.е. на его вход приходит сразу три вершины, принадлежащие одному треугольнику. Так же этот шейдер способен генерировать новые произвольные (по заданному пользователем алгоритму) вершины и удалять старые. На этом этапе необходимо сначала проверить вершины на видимость камерой, и если необходимо пересчитать некоторые их них. Видимой для камеры вершиной в данном случае считается вершина с неотрицательной координатой z (вся полусфера). В случае если в треугольнике нет неотрицательных вершин, он не отправляется на дальнейшую отрисовку (вершины удаляются).

На этапе растеризации (рисунок 2) графический процессор закрашивает все пиксели, принадлежащие треугольнику, но при стандартном использовании графического конвейера применяется матрица проекций, трансформирующая координаты примитивов в проекционные. Координата z(глубина изображения) приводится к диапазону (0.0,1.0), а координаты x и y приводятся к виду, в котором видимая часть изображения имеет координаты (-0.5,0.5).

В сферической системе координат по оси x откладываем , по оси y - , r – выступает в роли глубины изображения, где - зенитный угол, - азимутный угол, r - кратчайшее расстояние до начала координат (камеры). Для компенсации отсутствия проекционной матрицы необходимо разделить r на некоторую, заранее большую максимального расстояния в исследуемой модели величину (таким образом, приведя к нужному диапазону); в 4 компоненту координаты вершины в сферической системе координат вписать ( так как диапазон сферических координат (-,)

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

Полученное изображение обрисовывается не на экран, а в текстуру, и подсчитывается количества пикселей принадлежащих каждой попавшей в кадр поверхности, и делению этой величины на полное число пикселей в круговой части изображения. Полученные числа аддитивно вносятся в таблицу угловых коэффициентов излучения для дальней обработки[7-9].

Для подсчёта угловых коэффициентов на основе изображения выгодно так же использовать ресурсы графического процессора, чтобы, во-первых, не тратить время на выгрузку каждой картинки из видеопамяти в память компьютера (разрешение изображения оптимально делать не менее 2000 на 2000 пикселе), во-вторых, графическое ядро справится с задачей обработки изображения значительно быстрей. Для решения поставленной задачи предлагается использовать новшество DirectX 11 – вычислительный шейдер (Compute Shader). Вычислительный шейдер это программа, заменяющая собой весь графический конвейер, позволяя использовать графическое ядро для общих задач.

3. Исследование погрешностей

Для исследования погрешностей, вносимых разработанным методом, рассмотрим простейший пример – полый куб, каждая из граней которого состоит из двух треугольников и является излучающей поверхностью. Исключим из расчетов излучение внешние поверхностей куба. В этом случае шесть внутренних граней куба представляют собой замкнутую систему, и при равных начальных значениях температур не должно быть перераспределение энергии (значения температур не изменяются во времени).

В модели[7-9], разработанной на основе зонального метода, применяется рекурсивный учет отражения и одним из критериев корректности матрицы излучения является нулевая сумма элементов каждой строки (иначе в замкнутой системе при равно температуре всех её элементов будет происходить перераспределение температур, чего быть не должно). Назовем дефицитом излучения модели отношение отклонения суммы элементов строки от нуля к величине диагонального элемента. По сути дефицит излучения отражает количество излучения «потерянного» излучающей ячейкой из-за принятых в модели допущений. Ниже представлены таблицы коэффициентов дефицита излучения (в процентах) для куба при различных уровнях тесселяции и глубине отражения (в случае куба все 6 граней равнозначны).

Таблица 1: Коэффициенты дефицита излучения при ε = 0.2

Глубина отражения	1	2	3	4	5	10	20	30	50	100
Тесселяция
1	85.00	83.14	82.32	81.96	81.81	81.68	81.68	81.68	81.68	81.68
3	72.53	64.87	59.51	55.72	53.07	47.97	47.02	47.00	47.00	47.00
5	69.05	59.01	51.31	45.33	40.75	29.85	26.47	26.27	26.26	26.26
10	67.14	55.67	46.44	38.94	32.91	16.82	10.32	9.73	9.67	9.67
20	66.50	54.53	44.74	36.67	30.07	11.79	3.72	2.89	2.79	2.79
40	66.32	54.21	44.27	36.03	29.28	10.35	1.79	0.87	0.76	0.76
60	66.29	54.15	44.18	35.90	29.12	10.06	1.40	0.46	0.35	0.35

Таблица 2: Коэффициенты дефицита излучения при ε = 0.8

Глубина отражения	1	2	3	4	5	10	20	30	50	100
Тесселяция
1	51.70	51.20	51.15	51.14	51.14	51.14	51.14	51.14	51.14	51.14
3	19.28	17.20	16.84	16.78	16.77	16.76	16.76	16.76	16.76	16.76
5	10.79	8.07	7.55	7.45	7.44	7.43	7.43	7.43	7.43	7.43
10	6.23	3.11	2.50	2.38	2.35	2.35	2.35	2.35	2.35	2.35
20	4.71	1.45	0.80	0.67	0.65	0.64	0.64	0.64	0.64	0.64
40	4.29	0.10	0.34	0.20	0.18	0.17	0.17	0.17	0.17	0.17
60	4.21	0.91	0.25	0.11	0.08	0.08	0.08	0.08	0.08	0.08

Как видно из таблиц 1 и 2 при тесселяци больше 10 основную погрешность вносит недостаточное количество отражений, особенно при низком коэффициенте поглощения (при последнем вызове рекурсивной функции отраженная часть излучения теряется). Если при последнем вызове не учитывать отражение, то коэффициент дефицита резко уменьшится.

Вывод

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

Литература

1. Невский А.С. Лучистый теплообмен в печах и топках. 2-е изд., испр. и доп. М.: Металлургия,1971.440 с.

2. Спэрроу Э.М., Сесс Р.Д. Теплообмен излучением. Л.: Энергия, 1971. 296 с.

3. Зигель Р., Хауэлл Дж. Теплообмен излучением. М.: Мир, 1975. 936 с.

4. Юдаев Б.Н. Теплопередача. - Москва: Высшая школа,1973. - 360с.

5. Исаченко В.П., Осипова В.А., Сукомел А.С. Теплопередача.- Москва: Энергия,1975. - 480с.

6. Суринов Ю.А. Об итерационно зональном методе исследования и расчета локальных характеристик лучистого теплообмена. - Теплофизика высоких температур, 1972.Т.10. № 4. С. 844 – 852.

7. Лебедев, И.А. Модель термообработки металлических изделий в вакууме / И.А. Лебедев, А.П. Пономарев // VII Международная научно-практическая конференция «Энергосберегающие технологии в промышленности. Печные агрегаты. Экология»: сб. материалов. – Москва: МИСиС, 2014. – С. 270-280.

8. Лебедев, И.А. Управление равномерностью нагрева изделий в вакуумных электрических печных агрегатах / И.А. Лебедев, Д.В.Шатов // XII Всероссийское совещание по проблемам управления: сб. трудов. – Москва: ИПУ РАН, 2014. – С. 4315-4324.

9. Лебедев, И.А. Компьютерная модель нагрева теплозащитного контейнера в вакуумных печах / И.А. Лебедев, А.П. Пономарев // VI Международная научно-практическая конференция «Печные агрегаты и энергосберегающие технологии в промышленности. Экология»: сб. трудов. – Москва: МИСиС, 2012. – С. 319-320.