Моделирование теплообмена излучением: угловые коэффициенты излучения

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

Это третья статья из серии публикаций о моделировании теплообмена излучением в прозрачной среде. Прочтите первую и вторую части на нашем сайте.

Угловые коэффициенты излучения: краткое введение

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

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

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

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

Рассмотрим первый объект и оценим полный поток излучения, который он испускает. Часть испускаемого излучения уйдёт в окружающую среду, а часть попадёт на второй объект. Теперь самое время ввести понятие углового коэффициента излучения. Эта величина, которую мы обозначим F_{12}, определяется как доля излучения с поверхности 1 \left( A_{1} \right), которое попало на поверхность 2 \left( A_{2} \right). При условии однородного диффузного отражения и отсутствия затеняющих объектов и поверхностей угловой коэффициент излучения между поверхностями можно определить так:

Если в системе больше двух поверхностей, вполне возможно, что все они в зоне прямой видимости, поэтому угловой коэффициент F_{i j} будет соответствовать каждой паре поверхностей, где i,j — это индексы для каждой из N поверхностей, обменивающихся теплотой. Для любой пары поверхностей справедливо условие взаимности: A_i F_{i j} = A_j F_{ji}.

Заметим, что если поверхность вогнутая, то F_{i i} >0. Также стоит отметить, что плотность теплового потока в окружающую среду определяется угловым коэффициентом F_{i \rightarrow amb} = \left( 1 – \sum_{j=1}^N F_{ij}\right). Для замкнутых полостей угловой коэффициент излучения в окружающую среду равен нулю.

Три метода расчёта теплообмена излучением

В среде COMSOL Multiphysics реализовано три метода расчёта теплообмена излучением:

1. Метод прямого интегрирования
2. Метод полукуба
3. Метод трассировки лучей

1. Метод прямого интегрирования

Суть метода прямого интегрирования состоит в расчёте двойного интеграла по всем парам поверхностей. Метод применим только если между поверхностями нет препятствий или затеняющих объектов. Точность этого метода определяется исключительно порядком интегрирования.

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

2. Метод полукуба

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

В рамках метода полукуба для расчёта потока излучения внешние поверхности проецируются на набор разбитых на пиксели граней.

Затенение площадок эффективно учитывается с помощью z-буферизации, поэтому вычислительные затраты невелики. Единственным параметром этого метода является разрешение Radiation Resolution, которое определяет число пикселей. Точность выполнения условия взаимности повышается при увеличении разрешения, а угловой коэффициент излучения в окружающую среду всегда будет равен единице.

3. Трассировка лучей

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

Иллюстрация метода трассировки лучей в терминах дискретизации 3D полусферы при параметре разрешения, равном 4. Все 16 секций, расположенные в основании в шахматном порядке (слева), имеют одинаковую площадь. Стрелки направлены по направлению к вершине каждой секции на полусфере (справа).

При настройке метода трассировки лучей помимо порядка элементов используется ещё шесть параметров. Самым важным из них является разрешение Radiation Resolution, которое определяет начальное распределение лучей на полусфере (в 3D) или полукруге (в 2D), как показано на рисунке выше, где разрешение равно n_{res}=4.

Реализация метода начинается с разбиения внешнего пространства на n_{res}^2 секций в 3D (или на n_{res} в 2D) с последующим построением лучей, направленных в вершины каждой секции. Все эти секции имеют одинаковый угловой коэффициент, что означает, согласно аналогии Нуссельта, равенство площадей проекций этих секций на горизонтальную плоскость. В двумерных задачах, как показано ни рисунке ниже, полукруг разбивается на n_{res} секторов, и все они при проецировании на горизонталь формируют отрезки одинаковой длины. Отметим, что такое разбиение приводит к неоднородному по углу распределению лучей.

Иллюстрация метода трассировки лучей в 2D задаче. Каждый сектор полукруга (слева) при проекции на линии формирует отрезки одинаковой длины. Стреки направлены к границам каждого сектора (справа).

Трассировка луча позволяет определить плотность теплового потока в направлении луча и сравнить её с плотностью потока от соседних лучей. Если потоки, соответствующие соседним лучам, отличаются в пределах заданной точности, определяемой параметром Tolerance, тогда добавляются новые лучи, и процедура повторяется, пока не будет достигнуто максимальное число адаптаций, заданное параметром Maximum number of adaptations. Когда луч попадает на зеркальную или полупрозрачную поверхность, с этой поверхности испускается ещё один дополнительный луч. Максимальное число таких отражённых лучей задается параметром Maximal number of reflections. По умолчанию значение параметра равно 1000, что вполне достаточно, если только не моделируется множество отражений в полости с коэффициентом зеркального отражения, превышающим 0.99.

Параметр Angular dependent properties применяется в случае, когда в модели имеются поверхности с угловой зависимостью степени черноты. Опция по умолчанию Full для разрешения позволяет получить наиболее точные результаты, но при этом её использование ведёт к повышенным затратам вычислительных ресурсов по сравнению с опцией Interpolation function, которая позволяет задать, насколько точно нужно аппроксимировать угловую зависимость.

Чтобы получить достоверные результаты моделирования, прежде всего необходимо проанализировать параметры Radiation resolution и Maximum number of adaptations, поэтому важно понимать взаимосвязь между этими настройками. Давайте рассмотрим двумерный пример, в котором лучи испускаются элементом, расположенным в центре. Следует отметить, что этот график приведен только для визуализации; расчётные лучи сами по себе не отображаются. Мы рассмотрим внешнюю полусферу со степенью черноты, равной единице (что эквивалентно нулевому коэффициенту отражения) при фиксированной однородной температуре. В этом случае каждый луч соответствует одинаковой плотности излучения, и даже при минимальном разрешении мы сможем получить корректные результаты для плотности потока. Повышение разрешения (увеличение числа лучей) не приведёт к увеличению точности, а процедура адаптации не будет задействована.

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

Добавление малого объекта, который лежит на одном из лучей, приведет к запуску процедуры адаптации лучей в близлежащем пространстве. Чем выше число адаптаций — тем выше точность.

Наконец, добавим еще один малый объект в точке, которая не лежит ни на одном из начальных лучей. В этом случае число максимальных адаптаций будет не важно. Этот второй объект не станет “виден” ни одному из исходных лучей. Чтобы заметить этот второй объект, необходимо также повысить разрешение.

Если малый объект не “виден” ни одному из исходных лучей вследствие низкого разрешения, процедура адаптации не запускается, и объект исчезнет. В этом случае необходимо повысить разрешение.

Использование групп излучения

Дополнительно ко всем перечисленным выше методам можно использовать так называемые группы излучения. Можно существенно снизить вычислительные затраты, особенно в моделях с несколькими отдельными полостями, если задать группы поверхностей, которые могут “светить” только друг на друга. Группы следует использовать очень внимательно, поскольку в случае неверного выбора можно получить некорректные результаты.

Если некоторые группы поверхностей в принципе не могут “светить” друг на друга, целесообразно воспользоваться функционалом группировки Groups. Слева разными цветами показаны разумно сформированные группы. Справа пример менее разумного применения групп.

Другие настройки

В моделях с подвижными объектами необходимо пересчитывать угловые коэффициенты. Частота обновления значений угловых коэффициентов задаётся параметром View factors update threshold. По умолчанию значения обновляются на каждой итерации нелинейного решателя. Этот способ даёт наиболее точный результат, но при значительных вычислительных затратах. В случае, когда в результате перемещения или деформации объекта угловые коэффициенты изменяются незначительно, можно полностью отключить процедуру обновления их значений. Также можно настроить периодическое обновление значений, задав временной интервал или условное выражение.

Настройка процедуры расчёта угловых коэффициентов в интерфейсе Surface-to-Surface Radiation.

Кроме того, в разделе View Factor можно настроить сохранение значений угловых коэффициентов на жёсткий диск. Эта опция позволит ускорить решение масштабных моделей, однако размер файла модели на диске может значительно увеличиться, особенно при использовании метода полукуба. Использовать эту опцию целесообразно только в случае, когда геометрическая модель не изменяется.

Параметр Geometry Representation играет роль, если в модели для описания геометрии используются нелинейные функции формы. Если повысить порядок дискретизации, то данная опция влияет на учёт кривизны элементов.

Наконец, для параметра Jacobian Contribution по умолчанию задана опция Only local contributions to radiosity. Этот параметр, появившийся в COMSOL Multiphysics® версии 6.0, позволяет кратно снизить требования к оперативной памяти и ускорить вычисления. Однако в случае чисто радиационного охлаждения и значительной разности температур поверхностей задача может не сойтись. Если наблюдается отсутствие сходимости, измените настройку параметра на Include contributions from total irradiation.

Расчёт и визуализация угловых коэффициентов

О том, как рассчитать взаимные угловые коэффициенты излучения, если в модели имеется несколько поверхностей, можно узнать из ранее опубликованной статьи блога. Также иногда может быть полезно визуализировать значения угловых коэффициентов излучения с одной поверхности на все остальные элементы модели. Это можно сделать с помощью оператора elemint(order,expression), который выполняет интегрирование по Гауссу по каждому элементу. Порядок интегрирования задается первым аргументом оператора, а чтобы вычислить угловой коэффициент, нужно воспользоваться операторами radopu() и radopd() в качестве второго аргумента. Например, выражение

elemint(1,comp1.rad.radopu(S1,0))/intS1(1)/dvol

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

Переменная S1 должна быть задана равной единице на тех поверхностях, которые “светят”, и нулю на тех, которые “освещаются”. Ниже показан пример визуализации, построенной с помощью такого выражения. Дополнительно поделив размер освещенных элементов на переменную dvol, мы получим график, который показывает интенсивность освещенности с заданного набора поверхностей, как показано ниже.

Визуализация освещенности двух частично затеняющих друг друга блоков излучением со сферы.

Заключительные замечания

Мы рассмотрели три ключевых вопроса моделирования теплообмена излучением между поверхностями в прозрачной среде. Сначала мы обсудили, как поверхности испускают тепловое излучение. Затем мы показали, как падающее излучение поглощается, отражается и пропускается поверхностью. И теперь мы рассказали об угловых коэффициентах и методах их расчёта. Теперь мы с уверенностью может приступать к моделированию теплообмена излучением!