Практическая значимость

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

1. Зависимости теплофизических свойств материалов от температуры и времени.

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

3. Зеркальной компоненты отражения от экранных поверхностей.

4. Непланковского спектра излучения элементов конструкции.

5. Изменения положения внешних источников, связанного с движением обсерватории по орбите.

Материалы диссертации обсуждались на:

1. Международной конференции по оболочкам и пространственным конструкциям IASS-2007;

2. 8-ой Международной конференции и выставке «Системы проектирования, технологической подготовки производства и управления этапами жизненного цикла промышленного продукта» CAD/CAM/PDM – 2008;

3. XII Школе молодых ученых «Актуальные проблемы физики» и II Школе-семинаре «Инновационные аспекты фундаментальных исследований»;

4. Шестнадцатой Международной конференции по Вычислительной механике и современным прикладным программным системам ВМСППС’2009;

5. II Всероссийской конференции «Актуальные проблемы ракетно-космического приборостроения и информационных технологий»,

6. Научно-техническом семинаре в НПО им. С.А. Лавочкина,

7. Объединенном семинаре ИМП им. М.В. Келдыша и МГТУ им. Н.Э. Баумана «Механика и управление в робототехнических системах».

8. Всероссийской астрономической конференции ВАК-2010.

9. V Международной конференции «Параллельные вычисления и задачи управления» (PACO’2010).

Основные положения диссертации опубликованы в следующих работах:

1. Шаенко А.Ю., Милютин Д. С. Теплообмен в радиационных экранах больших космических телескопов // Доклады Академии Наук. – 2010. Том. 431. – с. 621-624.

2. Шаенко А.Ю., Милютин Д. С. Нестационарный радиационно-кондуктивный расчет больших космических обсерваторий // Вестник компьютерных и информационных технологий. – 2010. №9. – с. 3-6.

3. Шаенко А.Ю. Методика расчета теплообмена в радиационных экранах больших космических обсерваторий // Космонавтика и ракетостроение. – 2011. №1(62). – с. 57-64.

4. Bujakas V., Dmitriev V., Shaenko A. Deployable radiation screen for large space telescope. Proc. of International Symposium on spatial structures IASS- 2007. – 2007. –P.1-8.

5. Шаенко А.Ю. Распределенный параллельный расчет радиационно-кондуктивного теплообмена методом Монте-Карло на базе графических ускорителей // Труды V Международной конференции «Параллельные вычисления и задачи управления» (PACO’2010). – 2010. – с. 281-294.

6. Шаенко А.Ю. Основные результаты теплового расчета космической обсерватории "Миллиметрон" новым методом // Труды Всероссийской астрономической конференции ВАК-2010. – 2010. – с. 26.

7. Шаенко А.Ю. Разработка новой методики теплового расчета радиационных экранов больших космических телескопов // Труды Шестнадцатой Международной конференции по Вычислительной механике и современным прикладным программным системам ВМСППС’2009. – 2009 – с. 757.

8. Шаенко А.Ю. Разработка методики расчета лучистого теплообмена в радиационных экранах космической обсерватории // Труды 8-ой Международной конференции и выставки «Системы проектирования, технологической подготовки производства и управления этапами жизненного цикла промышленного продукта» CAD/CAM/PDM – 2008. – 2008 – с. 168-170.

9. Буякас В.И., Троицкий В.Ф., Шаенко А.Ю., Гордиенко А.М., Гришин Н.С. Моделирование задач укладки и раскрытия радиационных экранов большого космического телескопа // Труды 8-ой Международной конференции и выставки «Системы проектирования, технологической подготовки производства и управления этапами жизненного цикла промышленного продукта» CAD/CAM/PDM – 2008. – 2008 – с.164-167.

10. Шаенко А.Ю. Разработка методики расчета лучистого теплообмена в радиационных экранах космической обсерватории // Тезисы XII Школы молодых ученых «Актуальные проблемы физики» и II Школе-семинаре «Инновационные аспекты фундаментальных исследований». – 2008 – с. 73-74.

1. Тепловой расчет радиационных экранов больших космических телескопов (ОБЗОР)

Рост потребности астрономов в обсерваториях с высоким угловым разрешением и высокой чувствительностью приемной аппаратуры приводит к необходимости разрабатывать как телескопы с все большей апертурой, так и средства повышения чувствительности инструментов, в частности, обеспечивать их охлаждение [7, 23, 79, 86, 87]. Использовать уже отработанные технологии охлаждения на основе испарения жидкого гелия для крупногабаритных космических конструкций не представляется возможным, так как для этого потребуются чрезмерно большие запасы хладагента. Поэтому на перспективных космических обсерваториях с большими размерами зеркал все шире применяются радиационные экраны – устройства, обеспечивающие альтернативный способ охлаждения и поддержания теплового режима. Применение радиационных экранов обычно подразумевает ограничения на режимы ориентации космических аппаратов с целью постоянного затенения защищаемого объекта экраном.

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

Рис. 1.1. Экранно-вакуумная теплоизоляция. Слои металлизированной пленки Mylar проложены тканью типа «газ»

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

1.1. Радиационные экраны

Исторически первой космической обсерваторией с радиационным экраном можно считать Astronomical Netherlands Satellite (ANS) [30, 31], оснащенный аппаратурой для съемки в мягком и жестком рентгеновских и ультрафиолетовом диапазонах. Эта обсерватория была выведена на околоземную солнечно-синхронную орбиту в 1974 году. Ее радиационный экран закрывал часть корпуса, постоянно обращенную к Солнцу, при этом поддерживалась постоянная ориентация экраном на Солнце. Экран был выполнен нераскрывающимся. Внешний вид обсерватории ANS показан на Рис. 1.2. Схема аппарата изображена на Рис. 1.3.

Рис. 1.2. Внешний вид обсерватории ANS	Рис. 1.3. Схема обсерватории ANS

Радиационный экран аппарата Helios-1 [28], выведенного на околосолнечную эллиптическую орбиту в 1974 году, был предназначен для уменьшения тепловой нагрузки на внутренние системы во время близких пролетов Солнца. Первый пролет Helios-1 совершил 15 марта 1975 года, пройдя на расстоянии 47 млн. км от центра Солнца. Во время пролета аппарат вращался относительно продольной оси. Особенностью радиационного экрана Helios-1 является использование жестких металлических зеркал в качестве отражателей, в отличие от металлизированных пленок на других обсерваториях. Внешний вид Helios-1показан на Рис. 1.4.

Рис. 1.4. Внешний вид Helios-1

В дальнейшем радиационные экраны применялись в основном на космических аппаратах, предназначенных для работы в радиодиапазоне и для изучения анизотропии реликтового излучения, таких как COBE (запуск в 1989 году) [56], Odin (запуск в 2001 году) [21], WMAP (запуск в 2001 году) [9], Plank (запуск в 2009 году) [19]. Так же экраны применялись на обсерваториях инфракрасного диапазона с большими размерами зеркала, например, на Космическом телескопе Спитцера (запуск в 2003 году) [20], ASTRO-F (запуск в 2006 году) [35], Космической обсерватории Гершеля (запуск в 2009 году) [33]. Типичная инфракрасная обсерватория с нераскрывающимся радиационным экраном, Космический телескоп Спитцера, показана на Рис. 1.5. Схема обсерватории показана на Рис. 1.6. Типичный космический аппарат для исследования анизотропии реликтового излучения, WMAP, показан на Рис. 1.7. Схема внутреннего устройства WMAP изображена на Рис. 1.8.

Рис. 1.5. Внешний вид Космического телескопа Спитцера	Рис. 1.6. Схема Космического телескопа Спитцера
Рис. 1.7. Внешний вид космического аппарата WMAP	Рис. 1.8. Схема космического аппарата WMAP

Особенностью конструкции инфракрасных обсерваторий является нераскрывающийся экран, в некоторых случаях объединяемый с солнечными батареями. Радиационные экраны космических аппаратов, работающих в радиодиапазоне, отличаются от экранов инфракрасных обсерваторий в первую очередь раскрываемостью и назначением. Если экраны на Космическом телескопе Спитцера, аппарате ASTRO-F и Космической обсерватории Гершеля предназначены в основном для охлаждения приемной аппаратуры, то на аппаратах COBE, Odin, WMAP и Plank они дополнительно играют роль бленд, ограждая приемники от помех.

Достигнутый уровень развития технологии радиационных экранов позволяет поддерживать температуру защищаемых объектов на уровне 60К при использовании только экранов при соблюдении необходимой ориентации. Более низкая температура, вплоть до уровня 0,1К, поддерживается в течение срока продолжительностью до 15 месяцев с помощью совместной работы холодильных машин и экранов. Указанные параметры относятся к микроволновой обсерватории Plank, обладающей наиболее совершенной системой охлаждения в настоящее время.

Требования к перспективным обсерваториям более высокие. Экран обсерватории JWST должен будет обеспечивать температуру 40К зеркала диаметром 6,5 метров в течение 5,5 лет с возможностью продления до 10,5 лет [23]. При этом температуру инструментов необходимо будет поддерживать на уровне 7К с помощью холодильных машин. Внешний вид обсерватории JWST показан на Рис. 1.9, схема устройства – на Рис. 1.10.

Рис. 1.9. Внешний вид обсерватории JWST	Рис. 1.10. Схема обсерватории JWST

Обсерватория «Миллиметрон» будет работать в течение 3 лет с температурой главного зеркала 4К и температурой инструментов до 0,1К и затем еще 7-10 лет с общей температурой 50К [79]. Поддержание теплового режима зеркала и инструментов будет осуществляется комбинированной системой из холодильных машин и экранов. Охлаждение до 50К будет обеспечиваться только за счет экранов. Диаметр главного зеркала обсерватории – 12 метров. Схема обсерватории «Миллиметрон» показана на Рис. 1.11.

Рис. 1.11. Схема обсерватории «Миллиметрон»

Особенностями радиационных экранов перспективных космических обсерваторий являются:

- раскрываемость,

- большие линейные размеры,

- малая жесткость конструкции,

- возможность существенного искажения формы под действием нагрузки,

- работа при криогенных температурах.

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

Известно, что спектр излучения поверхностей реальных материалов может значительно отличаться от планковского [59], см. Рис.1.12. Расчет теплообмена между несерыми поверхностями представляет собой серьезную математическую проблему [93].

Рис. 1.12. Зависимость степени черноты поверхности, покрытой краской Cat-a-Lac, от длины волны [59]

Дополнительные сложности при проведении расчета создаются большой разницей температур экранов и аппаратуры обсерватории в транспортном и в рабочем состоянии. Например, для обсерватории «Миллиметрон» при температуре выведения, составляющей 300К, максимум спектра теплового излучения находится в районе 10 мкм, в то время как при 4К, рабочей температуре главного зеркала, максимум находится в районе 725 мкм. Термооптические коэффициенты поверхностей обсерватории для указанных диапазонов волн значительно отличаются [16]. Кроме того, термооптические и теплофизические характеристики материалов и покрытий существенным образом зависят от температуры [59, 10].

1.2. Методы расчета радиационного теплообмена

Расчет радиационного теплообмена основан на составлении уравнения теплового баланса для каждого излучающего и поглощающего тела. Для этого необходимо рассчитывать падающие на поверхность тела и излучаемые поверхностью потоки излучения. Способы их расчета зависят от набора предположений о свойствах поверхностей. Наиболее широко в литературе представлен радиационный расчет в предположении, что спектр излучения тел определяется законом Планка, и что поверхности тел отражают излучение диффузно [41, 84, 82, 89, 98, 93, 17, 47, 94, 102, 24, 78, 44, 70, 52, 85, 43, 91, 90, 97, 87, 86, 96]. В этом виде расчета вводится понятие угловых коэффициентов – поглощающихся в телах долей исходящих потоков излучения.

Следует отметить, что единое понятие «угловые коэффициенты» подразделяется на ряд более частных понятий в зависимости от принятых допущений о свойствах поверхностей. В случае расчета теплообмена между серыми диффузными поверхностями поглощенные доли энергии называются в русскоязычной литературе угловыми коэффициентами [85], в англоязычной – view factors [75]. Если при расчете учитывается зеркальное отражение, то в этом случае величина в англоязычной литературе имеет название – extended view factors [75], в русскоязычной литературе устоявшего обозначения данной величины нет. В наиболее сложном случае расчета – при учете и зеркального отражения и пропускания доли поглощенной энергии в русскоязычной литературе называются разрешающие угловые коэффициенты [85], в англоязычной – radiative exchange factors [75]. В настоящей работе все перечисленные понятия упоминаются как угловые коэффициенты с явным указанием принятых допущений.

Расчет значений угловых коэффициентов связан с двойным интегрированием по площадям выражений, зависящих от формы поверхностей, их взаимного расположения и ориентации [47, 93, 98, 5]. Получающиеся интегралы разрешаются в элементарных функциях только для некоторых сочетаний поверхностей, таких например, как «плоскость-плоскость», «двугранный угол», «вложенные сферы». Задача несколько облегчается наличием определенного количества известных решений для некоторых сочетаний поверхностей, но в большинстве случаев возникает необходимость проводить интегрирование численно.

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

- прямое численное интегрирование выражений для угловых коэффициентов,

- геометрические методы,

- методы Монте-Карло.

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

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

Рис. 1.13. Бесконечно малая площадка и построенная над ней полусфера [75]

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

Методы полусферы [75] основаны на аппроксимации полусферы над площадкой различными геометрическими фигурами: плоскостями, пирамидами, прямоугольными призмами [14] (так называемым полукубом, от английского hemicube), многогранниками, что позволяет упростить построение проекций. Важным моментом при использовании методов полусфер является применение пикселей (от английского pixel, picture element) – дискретов поверхности аппроксимирующей фигуры. Их применение позволяет упростить вычисление доли площади, занимаемой проекцией поверхности на полусфере, построенной вокруг площадки.

Метод Монте-Карло [51, 62] расчета угловых коэффициентов [37, 53, 64, 67, 93] основан на физической сущности процесса радиационного теплообмена. С физической точки зрения угловые коэффициенты между i -ой и j -той поверхностями являются отношением количества фотонов, поглощенных j -той поверхностью, к общему количеству фотонов, излученных i -ой поверхностью. В большинстве случаев число фотонов, излученных поверхностью за единицу времени, является чрезмерно большим для непосредственного моделирования хода каждого фотона. К тому же, в большинстве случаев точный расчет траектории отдельного фотона не возможен в принципе, что приводит к необходимости использовать статистические методы расчета хода излучения. Для этого весь набор фотонов, излученных с заданной поверхности за единицу времени, разбивается на заданное количество групп, пучков, и определяется доля мощности излучения, переносимая каждой группой. Затем для каждого пучка случайным образом выбирается направление излучения и определяется, попадает ли пучок в какую-либо поверхность и если попадает, то в какую именно. Для каждой поверхности ведется учет попавших пучков и определяется суммарная доля полной мощности излучения.

Методы полусферы и Монте-Карло широко используются для создания программных комплексов расчета радиационного теплообмена [2, 8, 1, 71, 26, 74, 27, 39, 73, 72, 34, 54, 69, 46, 3, 92]. Множество таких программ построено на допущениях о телах как серых излучателях и учитывают только диффузный характер отражения. В большинстве комплексов учтена зависимость теплопроводности и теплоемкости от температуры и/или времени. Однако в применении к расчету радиационных экранов, работающих как при комнатных, так и при криогенных температурах, указанные допущения приводят к чрезмерно большой погрешности, что вынуждает отказаться от допущения о планковском спектре излучения тел, учитывать зависимость термооптических свойств от температуры, времени и длины волны падающего излучения и принимать во внимание зеркальное отражение и пропускание излучения элементами конструкции.

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

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

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

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

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

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

1.3. Статистические методы расчета угловых коэффициентов

В случае диффузных серых поверхностей метод Монте-Карло расчета угловых коэффициентов реализуется следующей последовательностью действий, повторяемой для каждой поверхности [93]:

- на i -й поверхности случайным образом выбирается точка;

- из выбранной точки в случайном направлении излучается пучок фотонов;

- проверяется пересечение пучка с поверхностями исследуемого объекта;

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

- если пучок поглощается, счетчик пучков j -й поверхности увеличивается на единицу и происходит переход выбору точки и излучению следующего пучка;

- если пучок отражается диффузно, то отраженный пучок излучается в случайном направлении;

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

- цикл повторяется заданное количество раз для каждой поверхности;

- угловой коэффициент теплообмена между i -й и j -й поверхностями находится как отношение количества пучков, попавших в j -ю поверхность, к количеству пучков, излученному i -й поверхностью.

В случае учета более сложных видов взаимодействия, зеркального отражения и пропускания, и зависимости термооптических характеристик от длины волны и направления алгоритм необходимо записать иначе [93]:

- на i -й поверхности случайным образом выбирается точка;

- случайным образом в соответствии с заданным распределением определяется длина волны пучка;

- из выбранной точки в случайном направлении излучается пучок фотонов;

- проверяется пересечение пучка с поверхностями исследуемого объекта,

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

- если пучок поглощается, счетчик пучков j -й поверхности увеличивается на единицу и происходит переход к выбору точки и излучению следующего пучка;

- если пучок отражается зеркально, то в соответствии с законом зеркального отражения строится новое направление распространения излучения и происходит переход к излучению пучка;

- если пучок отражается диффузно, то случайным образом выбирается новое направление распространение излучения и происходит переход к излучению пучка;

- если пучок проходит сквозь поверхность, направление его распространения не изменяется;

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

- цикл повторяется заданное количество раз для каждой поверхности;

- угловой коэффициент теплообмена между i -й и j -й поверхностями находится как отношение количества пучков, попавших в j -ю поверхность, к количеству пучков, излученному i -й поверхностью.

Описанные реализации метода являются наиболее распространенными и, соответственно, встречаются во множестве вариантов. Основные модификации касаются двух пунктов приведенных алгоритмов – излучения пучка и проверки его пересечения с поверхностями объекта. Обе они направлены на повышение производительности расчета.

Изменения процесса излучения в основном направлены на повышение равномерности распределения пучков в исследуемом объекте. Среди этих модификаций можно выделить метод Маллея (Malley’s) [75, 48], метод косинусов [75] и метод стратифицированной полусферы (stratified hemisphere) [75].

Метод Маллея заключается в следующем:

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

- на круге случайным образом распределяются точки в количестве равном количеству излучаемых пучков;

- вокруг центра круга строится полусфера также единичного радиуса;

- точки с круга параллельно проецируются на поверхность полусферы;

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

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

Метод стратифицированной полусферы реализуется следующей последовательностью действий:

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

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

- на каждом дискрете случайным образом выбирается точка;

- вокруг центра круга строится полусфера единичного радиуса;

- точки с круга параллельно проецируются на поверхность полусферы;

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

Второй областью модификации исходного метода Монте-Карло является построение пересечения пучка с поверхностью.

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

Задача определения пересечения луча с поверхностью произвольной формы решается не только при расчете радиационного теплообмена, но и при построении реалистичных трехмерных изображений [5, 103, 75, 67, 14, 12, 22]. Более того, при подробном рассмотрении эти области можно назвать родственными.

Задача синтеза реалистичных трехмерных изображений обычно формулируется следующим образом [101]:

- заданы объекты, подлежащие построению (визуализации);

- задано расположение источников света;

- задано положение наблюдателя (камера) и направление его взгляда;

- на экране монитора необходимо построить изображение объектов такими, какими их видит наблюдатель.

Фактически, описанная задача сводится к определению цветов пикселей на экране монитора. Для ее решения предложено множество разнообразных алгоритмов, таких как z-буферизация [12, 103], метод многократных отражений [77, 15, 13], метод кажущейся светимости [5, 29, 37] и множество других [22, 25, 76]. В ряде таких методов решается задача расчета хода излучения, для чего привлекается математический аппарат, использующий угловые коэффициенты. Поэтому в области построения реалистичных трехмерных изображений большое количество работ также посвящено расчету численных значений угловых коэффициентов.

В одном из используемых методов, а именно в методе многократных отражений, в англоязычной литературе получившем название ray tracing [77, 13, 4], значения угловых коэффициентов определяются методом Монте-Карло, аналогичным методу, описанному в разделе 1.3. При реализации этой последовательности действий, так же как в расчете радиационного теплообмена, необходимо находить пересечения лучей с поверхностями объектов.

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

За прошедшее время в этой области достигнут существенный прогресс. Производительность современных видеоускорителей выросла более чем на порядок по сравнению с первыми образцами [57, 83]. Более того, сложилась парадоксальная ситуация – в некоторых задачах, связанных, например, с решением больших систем линейных алгебраических уравнений, видеоускорители оказались значительно более эффективны, чем центральные процессоры. Производители видеоускорителей, осознав сложившуюся ситуацию, приступили к выпуску специальных версии своих изделий, предназначенных исключительно для проведения математических расчетов, и разработке программного обеспечения, позволяющего приспособить для этого уже существующие видеоускорители. Ряд работ посвящен исследованию появившихся возможностей проведения ресурсоемких вычислений на ПЭВМ [83, 32, 63].

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

1.5. Программные комплексы расчета радиационного теплообмена

Одними из первых программных комплексов расчета радиационного теплообмена являются комплекс SINDA/G (Systems Improved Numerical Differencing Analyzer / Gaski) [8] и TRASYS (Thermal Radiation Analyzer System) [2, 26, 74], разработанные в середине 1960-х годов и начале 1970-х соответственно по заказу NASA. Программа SINDA/G предназначена для расчета теплообмена в общей постановке методом тепловых балансов [81]. Идея метода тепловых балансов заключается в записи теплового потока между элементами либо как функции разницы их температур в случае кондуктивного теплообмена, либо как функции разницы четвертых степеней их температур в случае радиационного теплообмена. Расчет угловых коэффициентов, определяющих тепловые связи в случае лучистого теплообмена, и поглощающихся потоков излучения от внешних источников производится в программе TRASYS.

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

Альтернативным способом расчета угловых коэффициентов для использования в SINDA/G является программный комплекс NEVADA (Net Energy Verification And Determination Analyzer) [54]. Комплекс позволяет учитывать пропускание, диффузное и зеркальное отражение, задавать двунаправленную функцию отражения. NEVADA не позволяет задавать зависимость термооптических коэффициентов от длины волны. Расчет значений угловых коэффициентов производится методом многократных отражений, комбинированным с методом Монте-Карло.

Наиболее современным программным комплексом, используемым в настоящее время NASA для расчета угловых коэффициентов и потоков излучения, является TSS (Thermal Synthesizer System) [27, 71]. TSS реализует метод статистических многократных отражений, и позволяет рассчитывать зеркальное и диффузное отражение, и пропускание. Возможен учет переменности термооптических свойств.

Еще одним способом подготовки данных для расчета лучистого теплообмена в SINDA/G является THERMICA [69, 72]. Расчет угловых коэффициентов производится методом многократных отражений со статистическими испытаниями. При проведении расчета учитываются зеркальное и диффузные отражения, преломление и пропускание. Необходимо отметить, что THERMICA позволяет представить результаты расчета в виде, пригодном для использования в программном комплексах ESATAN [39] и THERMISOL [73]. Эти комплексы являются аналогами комплекса SINDA/G. ESATAN и THERMISOL были разработаны для нужд Европейского космического агентства.

В комплексе ESATAN заложен тот же принцип расчета, что и в SINDA/G – метод тепловых балансов [81]. Модулем для расчета радиационных связей между элементами для ESATAN является ESARAD [73, 34]. Для расчета угловых коэффициентов модулем используется метод многократных отражений со статистическими испытаниями. При расчете возможен учет диффузного и зеркального отражений и пропускания.

Фактически, ESATAN, SINDA/G и THERMISOL являются тепловыми решателями общего вида. Исходные данные для расчета радиационного теплообмена формируются для них в программных комплексах TRASYS, NEVADA, THERMICA и ESARAD.

Иной подход используется при расчетах в программных комплексах MSC.PATRAN/MSC.NASTRAN [46], ANSYS [3], ABAQUS [1] и NEi Thermal Advanced [36]. Перечисленные комплексы способны генерировать набор угловых коэффициентов и внешних тепловых потоков как самостоятельно, так и использовать ранее их значения, рассчитанные в TRASYS, THERMICA и аналогичных модулях, за исключением ABAQUS. Описанные комплексы позволяют производить расчет радиационного теплообмена с учетом и зеркального и диффузного отражений и пропускания как статистическим методом многократных отражений, так и методом полусфер.

Комплекс ABAQUS, не поддерживает возможность использования сторонних угловых коэффициентов. ABAQUS реализует расчет теплообмена исключительно в «диффузно-серой» постановке методом Монте-Карло с учетом переменности термооптических коэффициентов во время расчета.

Отечественным комплексом для анализа радиационного теплообмена, сертифицированным для расчета космических аппаратов, является пакет прикладных программ «ТЕРМ» [92] разработки ЦНИИМаш. Комплекс использует метод многократных отражений со статистическими испытаниями для расчета угловых коэффициентов и внешних тепловых потоков. При проведении расчетов угловых коэффициенты считаются неизменными, термооптические коэффициенты принимаются не зависящими от температуры.

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

Необходимо еще раз заметить, что упоминаемое в описании программных комплексов выражение «угловые коэффициенты» соответствует понятиям «разрешающие угловые коэффициенты» в русскоязычной литературе и radiative exchange factors в англоязычной.

1.6. Требования к методу расчета радиационного теплообмена, учитывающему сложные виды взаимодействия и переменные теплофизические и термооптические свойства

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

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

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

2. Расчет радиационно-кондуктивного теплообмена

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

2.1. Допущения и предположения

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

- решается нестационарная задача радиационно-кондуктивного теплообмена;

- исследуемая система представляется в виде набора трехмерных и двухмерных конечных элементов;

- материал трехмерных элементов не пропускает излучение;

- материал двумерных элементов достаточно тонкий и взаимодействует с излучением без преломления;

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

- распределение температуры и теплового потока по элементам на каждом шаге времени принимается равномерным;

- спектр излучения может определяться как законом Планка, так и задаваться явно;

- конвективный теплообмен отсутствует;

- материалы исследуемой системы считаются изотропными;

- теплоемкость и теплопроводность материалов зависят от температуры и времени;

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

- излучение взаимодействует с материалом посредством отражения (зеркального и диффузного), поглощения и пропускания;

- внешние источники излучения точечные, излучающие равномерно во все стороны;

- координаты внешних источников излучения зависят от времени;

- в качестве граничных условий задаются мощности тепловыделения в элементах и/или их температуры, зависящие от времени;

- излучение происходит равномерно со всей грани элемента;

- поглощение излучения происходит равномерно по всему объему элемента.

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

2.2. Математическая постановка задачи радиационно-кондуктивного теплообмена

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

,

(2.1)

где

- теплоприток в точку x за счет кондуктивного теплообмена; - теплоприток в точку x за счет лучистого теплообмена со всей областью V; - доля излучения, испущенного точкой y и поглощенного в точке x; , - мощность излучения, испускаемого из точек и соответственно; - мощность излучения, испускаемого источниками и поглощающегося в точке ; - мощность тепловыделения в точке от всех внутренних источников тепла; - мощность излучения, испускаемого m -ым источником и поглощающегося в точке x; - мощность внутреннего тепловыделения p -го внутреннего источника в точке x; M - количество источников излучения; P - количество внутренних источников тепла; V – исследуемая область; x, y - точки области; t - время; - температура в точках x, y соответственно, - удельная теплоемкость, - плотность, - теплопроводность.

Граничные условия уравнения (2.1) формулируются в виде

, , , ,

(2.2)

где

- заданная функция распределения температуры,

- область, в которой задается распределение температур,

- заданная функция распределения тепловыделения,

- область, в которой задается распределение тепловыделения.

Начальные условия уравнения (2.1) записываются в виде

,

(2.3)

где

- заданная функция распределения температуры,

- начальный момент времени,

.

Рис. 2.1 Исследуемая область

Для решения задачи теплового расчета, то есть для определения распределения температур в исследуемой области на заданном интервале времени, необходимо численно проинтегрировать уравнение (2.1) по времени с граничными условиями (2.2) и начальными условиями (2.3).

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

,

(2.4)

где - удельная теплоемкость i -го элемента; - его плотность; - объем элемента; - температура элемента; - время; - количество элементов, имеющих общую грань с i -м элементом; - мощность, передаваемая от k -го соседа i -го элемента посредством теплопроводности; - площадь i -го элемента; - плотность исходящего лучистого потока i -го элемента; - мощность, передаваемая от j -го элемента на i -й элемент посредством излучения; - количество элементов в модели; - мощность излучения, испускаемого m -ым внешним источником и поглощенного в i -ом элементе; - количество источников излучения; - тепловыделение от p -го внутреннего источника тепла в i -ом элементе, - количество внутренних источников тепла.

Задача нестационарного радиационно-кондуктивного теплового расчета состоит в численном интегрировании системы (2.4) при заданных начальных и граничных условиях на заданном интервале времени. В качестве начального условия выступает распределение температуры по конечным элементам в исходный момент времени (начальное «поле температур»). В качестве граничных условий – заданные температуры и тепловыделения в элементах конструкции и внешние тепловые потоки.

Ясно, что получить решение задачи (2.4) в аналитическом виде в общем случае невозможно. Для ее решения необходимо построить численный метод решения, включающий в себя математическую модель испускания и поглощения излучения, способ расчета кондуктивного теплообмена и метод интегрирования по времени. Так как на основе анализа литературы в качестве метода решения задачи был выбран метод Монте-Карло, то необходимо выбрать подходящий механизм получения случайных чисел и сформулировать математическую модель испускания и поглощения излучения в вероятностном виде, пригодном для проведения статистических испытаний.

2.3. Расчет испускаемых лучистых потоков

1.

2.

Тепловое излучение – электромагнитное излучение, испускаемое веществом, возникающее за счет внутренней энергии тела и имеющее непрерывный спектр. В простейшем случае, при использовании модели абсолютно черного тела (АЧТ), спектр теплового излучения описывается законом Планка [93]:

,

(2.5)

где

- мощность излучения с единицы площади поверхности в диапазоне частот излучения от до при температуре ,

- постоянная Планка,

- скорость света в вакууме,

- постоянная Больцмана.

Зависимость (2.5), представленная в графическом виде, имеет следующий вид (Рис. 2.1).

Рис. 2.1. Спектры излучения абсолютно черного тела при различных температурах

Поток теплового излучения c поверхности абсолютно черного тела на всех частотах рассчитывается в соответствии с законом Стефана-Больцмана:

,

(2.6)

где

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

.

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

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

.

(2.7)

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

.

(2.8)

Полученное выражение представляет собой кумулятивное распределение вероятности излучения фотона от частоты и температуры. Вид зависимости представлен на Рис.2.2.

Рис. 2.2. Кумулятивное распределение вероятности излучения фотона в зависимости от частоты и температуры АЧТ

Имея набор случайных чисел, равномерно распределенных на отрезке [0;1] и выражение (2.8), становится возможным получить распределение частот излучаемых фотонов, соответствующее спектру абсолютно черного тела при заданной температуре. Для этого случайные числа подставляют в левую часть выражения (2.8) и находят соответствующее им значение частоты.

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

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

- провести интегрирование спектра по всем частотам и получить аналог выражения (2.6):

,

(2.9)

- разделить на и проинтегрировать получившееся выражение по частоте:

,

(2.10)

Выражение (2.10) пригодно для построения частотного распределения фотонов излучения, соответствующего заданному спектру .

Помимо частоты фотоны также характеризуются направлением распространения. В простейшем случае, фотоны излучаются равномерно во все стороны. Равномерно распределенные угловые координаты получаются по методу Маллея [48] следующим образом:

- вокруг точки-центра излучения строится полусфера единичного радиуса;

- полусфера параллельно проецируется на излучающую поверхность;

- на проекции случайным образом распределяются точки, см. Рис. 2.3;

- точки с круга параллельно проецируются обратно на полусферу, см. Рис. 2.4;

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

Рис. 2.3. Точки, случайным образом распределенные по кругу

Рис. 2.4. Проекции точек на полусферу

В случае необходимости учета неравномерности излучения фотонов метод Маллея необходимо модифицировать. В частности, учет неравномерности может быть произведен изменением распределения точек по кругу. Выбирая неравномерное распределение точек по кругу, можно получить неравномерное распределение направления излучения фотонов, см. Рис. 2.5 и Рис. 2.6.

Рис. 2.5. Точки с неравномерным распределением по кругу

Рис. 2.6. Неравномерное распределение точек на полусфере

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

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

Рис. 2.7. Точки с неравномерным распределением по кругу для внешнего источника излучения

Рис. 2.8. Неравномерное распределение точек на полусфере для внешнего источника излучения

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

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

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

Плотности тепловых потоков и , исходящие с граней элементов, в выражении (2.4) равны величинам, рассчитываемым по формулам (2.6) и (2.9). Количество пучков, излучаемое элементами и источниками излучения на одном шаге времени задается в явном виде и обозначается и соответственно.

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

2.4. Расчет поглощаемых лучистых потоков

Излученный пучок фотонов, может либо выйти за пределы рассматриваемой области, либо поглотиться в элементе системы. Конечный результат распространения пучка определяется его первоначальным направлением распространения, его длиной волны, а также историей взаимодействия с элементами модели.

Распространение пучка в исследуемой системе связано с решением следующих подзадач:

- необходимо определить, существует ли точка пересечения пучка фотонов с элементами системы и если существует, то рассчитать ее координаты,

- в случае, если пересечение существует, то определить вид взаимодействия пучка с элементом,

- определить дальнейшее направление распространения пучка или провести учет его поглощения.

2.4.1. Расчет точки пересечения пучка с элементами модели

Координаты x, y и z точки пересечения пучка с элементом являются решением системы четырех алгебраических уравнений, состоящих из уравнений для пучка:

,

(2.11)

где