38712

МЕТОД ТЕПЛОВОГО РАСЧЕТА БОЛЬШИХ КОСМИЧЕСКИХ ТЕЛЕСКОПОВ И ЕГО ПРОГРАММНАЯ РЕАЛИЗАЦИЯ

Диссертация

Математика и математический анализ

ЛЕБЕДЕВА На правах рукописи Шаенко Александр Юрьевич МЕТОД ТЕПЛОВОГО РАСЧЕТА БОЛЬШИХ КОСМИЧЕСКИХ ТЕЛЕСКОПОВ И ЕГО ПРОГРАММНАЯ РЕАЛИЗАЦИЯ Специальность 05.3] МЕТОД ТЕПЛОВОГО РАСЧЕТА [0. Существующие программные комплексы теплового расчета предназначены в основном для расчета космических аппаратов с небольшим разбросом температур. Прямая реализация расчета космических аппаратов с большим разбросом температур и сложной геометрической конфигурацией по классической вычислительной схеме приводит к необходимости использовать суперЭВМ в то время как...

Русский

2013-09-29

4.57 MB

34 чел.

PAGE  2

РОССИЙСКАЯ АКАДЕМИЯ НАУК

ФИЗИЧЕСКИЙ ИНСТИТУТ ИМ. П.Н. ЛЕБЕДЕВА

На правах рукописи

Шаенко Александр Юрьевич

МЕТОД ТЕПЛОВОГО РАСЧЕТА

БОЛЬШИХ КОСМИЧЕСКИХ ТЕЛЕСКОПОВ

И ЕГО ПРОГРАММНАЯ РЕАЛИЗАЦИЯ

Специальность 05.13.18 – Математическое моделирование, численные методы и комплексы программ

Диссертация на соискание ученой степени

кандидата технических наук

Научный руководитель –

Доктор физико-математических наук

В.И. Буякас

Москва 2011

[0.0.0.1] РОССИЙСКАЯ АКАДЕМИЯ НАУК

[0.0.0.2] ФИЗИЧЕСКИЙ ИНСТИТУТ ИМ. П.Н. ЛЕБЕДЕВА

[0.0.0.3] МЕТОД ТЕПЛОВОГО РАСЧЕТА

[0.0.0.4] БОЛЬШИХ КОСМИЧЕСКИХ ТЕЛЕСКОПОВ

[0.0.0.5] И ЕГО ПРОГРАММНАЯ РЕАЛИЗАЦИЯ

[1] Заключение

[2]
Приложение. Краткое описание алгоритма T.H.O.R.I.U.M.

[3]
Список  литературы


Введение

Прогресс современной астрономии в настоящее время неразрывно связан с освоением новых диапазонов электромагнитного спектра, повышением чувствительности и углового разрешения инструментов. Решение этих задач стало возможным с выведением обсерваторий за пределы земной атмосферы. Так, например, телескоп SIGMA международной астрофизической обсерватории «Гранат» [], выведенной на орбиту в 1989 году, позволил впервые в мире произвести съемку области центра Галактики в жестком рентгеновском диапазоне. Построение современной космологической модели Вселенной стало возможным после запуска в 1990 году Космического телескопа Хаббла с беспрецедентным по тем временам угловым разрешением в 0,1 угловой секунды []. Первые в мире непосредственные спектроскопические наблюдения внесолнечных планет HD 189733b и HD 209458b были проведены при помощи инструмента MIPS чувствительностью в 1,5 мЯнских на инфракрасной обсерватории Спитцера [] , запущенной в 2003 году.

Одним из методов повышения чувствительности приемной аппаратуры является ее охлаждение. На действующих и завершивших работу обсерваториях, таких как IRAS [], ISO [], обсерваториях Спитцера и Гершеля [], охлаждение производится в основном с помощью криогенных систем. Конструктивно эти обсерватории представляют собой криостаты с размещенными внутри или инструментами, или телескопами целиком. Габаритные размеры зеркал таких обсерваторий ограничиваются диаметром головных обтекателей ракет-носителей, используемых для запуска обсерваторий на орбиту. Из перечисленных обсерваторий наибольшее зеркало диаметром 3,5 метра имеет обсерватория Гершеля, что близко к пределу возможностей современных средств выведения.

Разрабатываемые в настоящее время обсерватории JWST [], TPF-C [] и «Миллиметрон» [, , , ] будут иметь телескопы с главными зеркалами характерным размерами 6,5 метров, 8 метров  и 12 метров соответственно. Для размещения под головным обтекателем ракет-носителей эти телескопы необходимо будет выполнить раскрывающимися. Потребуется предусмотреть компактную укладку обсерватории в транспортном положении при размещении ее на носителе и обеспечить раскрытие обсерватории из транспортного положения в рабочее после выхода на орбиту.

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

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

С точки зрения конструкции радиационные экраны представляют собой тонкую металлизированную полимерную пленку типа лавсан (Mylar, Kapton), закрепленную на жестком каркасе. На обсерваториях ISO, IRAS, Спитцера и Гершеля каркасы экрана были выполнены неподвижными относительно корпуса телескопа. На разрабатываемых обсерваториях с раскрывающимися зеркалами, таких как JWST, TPF-C и «Миллиметрон», каркасы экранов потребуется также выполнять раскрывающимися.

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

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

Актуальность работы

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

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

Цель работы

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

Задачи исследования

Для достижения поставленной цели необходимо было решить ряд задач:

  1.  Разработать метод нестационарного радиационно-кондуктивного теплового расчета больших космических телескопов, учитывающий зеркальную составляющую отражения экранов, переменность теплофизических и термооптических характеристик материалов, непланковский спектра излучения элементов конструкции, а также переменный характер внешних источников, связанный с движением обсерватории по орбите;
  2.  Программно реализовать предложенный метод с учетом ограничений по  ресурсам ЭВМ;
  3.  Провести проверку достоверности и работоспособности программной реализации метода;
  4.  Используя разработанное программное обеспечение, оценить влияние точности изготовления, термооптических свойств материала геометрии экрана на температуру главного зеркала космического телескопа проекта «Миллиметрон».

На защиту выносятся:

  1.  Новый метод теплового расчета космических объектов сложной геометрии.
    1.  Программный модуль Т.H.О.R.I.U.M,  реализующий этот метод.
    2.  Результаты сравнительного анализа влияния геометрии, точности изготовления и теплофизических свойств материала экрана на температуру главного зеркала космического телескопа проекта «Миллиметрон», полученные с помощью Т.H.О.R.I.U.M.

Научная новизна работы заключается:

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

Достоверность результатов

Достоверность полученных результатов подтверждается поверочными расчетами тестовых примеров. В качестве тестовых примеров выбирались задачи с известными решениями. Хорошее совпадение результатов численных расчетов с известными решениями подтверждают достоверность метода.

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

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

  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.
  11.  
    Тепловой расчет радиационных экранов больших космических телескопов (ОБЗОР)

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

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

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

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

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

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

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

Рис. 1.3. Схема обсерватории ANS

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

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

В дальнейшем радиационные экраны применялись в основном на космических аппаратах, предназначенных для работы в радиодиапазоне и для изучения анизотропии реликтового излучения, таких как COBE (запуск в 1989  году) [], Odin (запуск в 2001 году) [], WMAP (запуск в 2001 году) [], Plank (запуск в 2009 году) []. Так же экраны применялись на обсерваториях инфракрасного диапазона с большими размерами зеркала, например, на Космическом телескопе Спитцера (запуск в 2003 году) [], ASTRO-F (запуск в 2006 году) [], Космической обсерватории Гершеля (запуск в 2009 году) []. Типичная инфракрасная обсерватория с нераскрывающимся радиационным экраном, Космический телескоп Спитцера, показана на Рис. 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 лет []. При этом температуру инструментов необходимо будет поддерживать на уровне 7К с помощью холодильных машин. Внешний вид обсерватории JWST показан на Рис. 1.9, схема устройства – на Рис. 1.10.

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

Рис. 1.10. Схема обсерватории JWST

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

                                                                              

,

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

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

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

,

(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) пригодно для построения частотного распределения фотонов излучения, соответствующего заданному спектру .

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

  •  вокруг точки-центра излучения строится полусфера единичного радиуса;
  •  полусфера параллельно проецируется на излучающую поверхность;
  •  на проекции случайным образом распределяются точки, см. Рис. 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). Количество пучков, излучаемое элементами и источниками излучения на одном шаге времени задается в явном виде и обозначается и  соответственно.

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

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

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

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

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

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

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

,

(2.11)

где

,, - координаты единичного вектора распространения пучка,

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

- параметр,

, ,  - координаты точки пересечения,

и уравнения плоскости элемента:

,

(2.12)

где

, ,  - координаты вектора единичной нормали к плоскости элемента,

, ,  - координаты точки на плоскости элемента,

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

  •  точка пересечения должна лежать внутри границ элемента,
  •  точка пересечения должна лежать по вектору распространения пучка, то есть .

Удобнее систему из уравнений (2.11) и (2.12) сначала разрешить относительно t, а затем найти координаты точки пересечения, подставив значение параметра в уравнение (2.11). Значение параметра t вычисляется по формуле:

.

(2.13)

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

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

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

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

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

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

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

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

Одним из наиболее широко распространенных картинных методов удаления невидимых граней, является метод z-буфера [, ]. Метод описывается следующей последовательностью действий:

  •  перед началом работы каждому пикселю присваивается значение глубины z = +∞;
  •  каждому пикселю ставится в соответствие глубина z как расстояние от картинной плоскости до грани;
  •  на картинную плоскость в растровой форме проецируются грани, затем для каждого пикселя грани находится свое значение глубины;
  •  на экран выводятся пиксели с наименьшей глубиной.

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

Наиболее распространенными программными средствами, позволяющими использовать ресурсы видеоускорителей, являются библиотеки OpenGL (Open Graphics Library) [, , , , ], которая разрабатывается некоммерческой организацией Khronos Group и Direct3D [] компании Microsoft.

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

OpenGL позволяет:

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

Direct3D решает задачи, аналогичные задачам OpenGL. Главной особенностью Direct3D является ее распространение исключительно в среде Microsoft Windows.

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

Основными объектами OpenGL, использованными в настоящей работе являются:

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

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

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

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

Этап расчета пересечения. Проводится на каждом шаге интегрирования и включает следующие операции:

  •  рассчитывается точка испускания пучка и направление его распространения согласно разделу 2.3 настоящей работы;
  •  в точку испускания пучка помещается узконаправленная камера, направление взгляда которой совпадает с направлением распространения пучка;
  •  в буфере изображения строится изображение модели, какой она выглядит из точки испускания пучка при взгляде в направлении распространения пучка;
  •  определяется цвет пикселя, находящегося в центре построенного изображения модели;
  •  в случае если цвет пикселя черный, то пучок считается покинувшим пределы рассматриваемой области системы;
  •  в случае если цвет пикселя отличается от черного, то пучок считается попадающим в элемент системы;
  •  в случае попадания пучка в элемент по таблице соответствия цвета номеру элемента определяется элемент, в который попадает пучок;
  •  по формулам (2.11)-(2.13) рассчитываются координаты точки пересечения пучка и элемента.

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

2.4.2.   Определение вида взаимодействия пучка с элементом

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

В общем случае направление распространения пучка после взаимодействия с элементом определяется с помощью двунаправленной функции рассеивания BSDF (Bidirectional scattering distribution function) [, ] – зависимости, связывающей доли энергии, отражаемые и пропускаемые в различных направлениях, с направлением падения излучения и его длиной волны. В самой общей постановке пучок может выходить из материала не обязательно в том же месте, где произошло взаимодействие, и преломляться без соблюдения закона Снелла. Примерный вид пучка до и после взаимодействия, описанного двунаправленной функцией рассеивания, показан на Рис. 2.9.

В более упрощенной форме взаимодействие пучка и элемента можно описать, разделив BSDF на двунаправленную функцию отражения BRDF (Bidirectional reflectance distribution function) и двунаправленную функцию пропускания BTDF (Bidirectional transmittance distribution function). В этом случае процессы преломления и пропускания происходят независимо. Примерный вид пучка до и после взаимодействия, описанного двунаправленной функцией отражения, показан на Рис. 2.10.

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

Рис. 2.9. Вид пучков  при использовании двунаправленной функции рассеивания BSDF

Рис. 2.10. Вид пучков  при использовании двунаправленной функции отражения BRDF

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

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

,

(2.14)

где

- значение функции в искомой точке ,

- значения функции в известных точках ,

- расстояние между  и ,

- заданный показатель степени.

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

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

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

,

(2.15)

где

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

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

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

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

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

Если величина не зависит от какого-либо параметра, то этот параметр исключается из вычисления расстояния между точками. Например, если величина не зависит от температуры, то тогда принимается .

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

  •  задать произвольный интервал для генерации случайных чисел, например интервал OD [0;1];
  •  разбить интервал OD на отрезки, пропорциональные величинам термооптических коэффициентов. Например, если для актуальных температуры, времени, длины волны и углов падения коэффициент зеркального отражения составляет, диффузного отражения , коэффициент пропускания равен и коэффициент поглощения , то полный интервал [0;1] необходимо разбить на следующие подынтервалы  [0;0,2),  [0,2;0,5),  [0,5;0,9) и  [0,9;1,0] соответственно, см. Рис. 2.11;
  •  определить конкретный вид взаимодействия путем генерации случайного числа на интервале OD. В случае если случайное число попадает в подынтервал , то считается, что пучок зеркально отражается, если число попадает в подынтервал , то пучок считается отраженным диффузно и т.д.

Рис. 2.11. Интервал OD и подынтервалы для определения вида взаимодействия

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

2.4.3.   Расчет направления распространения пучка после взаимодействия его с элементом

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

Рис. 2.12. Зеркальное отражение,  диффузное отражение и пропускание

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

2.4.4.   Учет поглощения пучка 

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

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

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

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

.

(2.16)

Соотношение (2.16) позволяет выразить в явном виде количество поглощенных пучков в зависимости от количества пропущенных

.

(2.17)

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

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

  1.   Безматричный распределенный расчет радиационного теплообмена

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

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

В литературе широко представлен способ хранения результатов расчета лучистых потоков в виде матрицы угловых коэффициентов F  [, , , , ]. Строка такой матрицы с индексом i представляют собой доли потока с i–го элемента, поглощенные во всей остальной модели. То есть, элемент матрицы  представляет собой часть потока излучения с i–го элемента, поглощенного в j–м элементе. В соответствии с принятыми обозначениями, значение элемента матрицы рассчитывается по формуле

,

(2.18)

где

- количество пучков фотонов, испущенных i–м элементом и поглощенных в j–м элементе,

- общее количество пучков фотонов, испущенных i –м элементом.

Значение мощности излучения  i–ого элемента, поглощающегося в j–м элементе, записывается следующим образом

,

(2.19)

где

- площадь l–ой грани i–ого элемента,

- количество граней i–ого элемента,

- поток излучения с l–ой грани в i–ого элемента.

Таким образом, суммарная мощность излучения, испущенного элементами модели и поглощенного в  j–м элементе на одном шаге времени, получается суммированием по i выражения (2.19)

,

(2.20)

где

- количество элементов в модели.

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

В случае же, если модель содержит большое количество элементов и/или угловые коэффициенты модели в процессе расчета меняются, то отдельно рассчитывать и хранить матрицу F становится не выгодно. Например, для модели, состоящей из 10000 элементов, объем оперативной памяти, необходимый для хранения матрицы угловых коэффициентов с точностью по 8 байт на один элемент матрицы, составляет 770 Мбайт. Для модели из 100000 элементов необходимый объем памяти составляет 76 Гбайт. Хранение матрицы такой размерности и операции с ней являются сложными задачами для современных персональных компьютеров.

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

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

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

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

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

.

(2.21)

Формула (2.21) получается из (2.20) путем подставки в нее выражения для расчета значения углового коэффициента (2.18). Заметим, что единственным членом в зависимости (2.21), определяемым посредством статистических испытаний, является  - количество пучков, поглощенных в j-м элементе.  Количество пучков, излучаемых на каждом шаге по времени , количество граней i–ого элемента  и их площади  задаются в качестве исходных данных. Потоки излучения с граней i–ого элемента  рассчитываются согласно выражениям (2.6) или (2.9).

То есть, статистические испытания необходимо проводить на каждом шаге по времени. Известно, что точность статистического расчета зависит от количества испытаний в степени 0,5. Отсюда следует, что, при сохранении точности расчета температурного поля в конце интервала интегрирования по времени, точность расчета выражения (2.21) может быть снижена в число раз,  пропорциональное количеству шагов по времени.

Запишем это утверждение в виде формулы. Стандартное отклонение случайной величины P, полученной по результатам  испытаний, вычисляется следующим образом:

 ,

(2.22)

где

- значение величины, полученное по результатам i–ого испытания,

- среднее значение величины.

Из выражения (2.22) следует, что увеличение допустимого значения стандартного отклонения на порядок позволяет на два порядка сократить количество испытаний.

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

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

Заметим, что мощность излучения от m-го источника излучения, поглощаемого в i-ом элементе можно определить по зависимости, сходной с зависимостью (2.21)

.

(2.23)

где

- количество пучков, излученных m-м источником и поглощенных l-ой гранью i-ого элемента,

 - общее количество пучков, излученных m-м источником,

 - суммарная мощность m-го источника излучения, задается в явном виде.

Используя выражения (2.21) и (2.23), можно вычислить мощность того излучения, которое испускается всеми элементами модели и источниками излучения, и поглощается в элементе, см. зависимость (2.4).

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

Например, если доступный при проведении расчета объем памяти ЭВМ составляет 1 Гбайт, то в нем можно хранить матрицу модели с примерно  элементами с точностью по 8 байт на элемент. В тоже время, при использовании безматричного метода количество элементов в модели может быть увеличено до .

2.5.2.   Распределенный расчет радиационного теплообмена

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

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

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

Сети персональных ЭВМ распространены значительно шире, при этом их суммарная вычислительная мощность может достигать весьма значительных величин. Так, например, один из наиболее известных проектов распределенных вычислений SETI@home сумел к 14 декабря 2009 года достичь вычислительной мощности 711,9 Терафлоп, что сравнимо с мощностью Cray XT5 (831,7 Терафлоп), занимающего третью позицию в рейтинге самых производительных суперкомпьютеров в мире.

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

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

  1.  Выделяется ЭВМ-сервер и задается ее адрес в локальной сети или Интернет.
    1.  На сервер устанавливается программное обеспечение в серверном варианте для проведения распределенных вычислений.
      1.  На ЭВМ-клиенты устанавливается программное обеспечение в клиентском варианте.

Этап проведения расчета. Проводится на каждом шаге интегрирования по времени и состоит в следующей последовательности действий:

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

Сервер производит предварительный анализ модели, приводит ее в вид, пригодный для проведения расчета.

Сервер формирует вектор температур на текущем шаге времени.

Клиенты через заданный промежуток времени проверяют наличие модели на сервере и при наличии модели производят ее считывание.

Клиенты производят оценку собственной производительности и сообщают о ее результатах серверу.

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

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

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

  1.   Расчет кондуктивного теплообмена

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

,

(2.24)

где

x - точка области, см. рис. 2.1;

t - время;

- температура в точке x;

- удельная теплоемкость в точке x в момент времени t при температуре T;

- плотность в точке x;

- теплопроводность в точке x в момент времени t при температуре T.

С учетом принятого допущения о равномерном распределении температуры по элементу правую часть уравнения (2.24) для i-го элемента можно переписать в следующем виде

,

(2.25)

где

- количество соседних элементов, имеющих общие грани с i-м элементом;

j-й элемент - j-й соседний элемент i-ого элемента, имеющий с ним общую грань;

 - средняя теплопроводность между центрами i-ого и j-ого элементов в момент времени t при температурах элементов  и  соответственно;

,  - расстояния от общей грани до центров i-ого и j-ого элементов соответственно;

 - расстояние между центрами i-ого и j-ого элементов;

 - площадь грани, общей для i-ого и j-ого элементов.

Средняя теплопроводность между центрами i-ого и j-ого элементов  рассчитывается по формуле

,

(2.26)

где

,  - теплопроводности i-ого и j-ого элементов.

Необходимо отметить, что расстояние  между центрами i-ого и j-ого элементов необходимо рассчитывать вдоль пути теплопередачи как сумму расстояний от центра одного элемента до общей грани и от общей грани до центра другого элемента, см. Рис. 2.13.

Рис. 2.13. Расстояние между центрами i-ого и  j-ого элементов

Фактически, выражения (2.25)-(2.26) реализуют метод тепловых балансов для расчета кондуктивного теплообмена [].

Используя выражения (2.25) и (2.26), можно рассчитать мощность, выделяющуюся в i-ом элементе за счет кондуктивного теплообмена с соседними элементами. Радиационная часть тепловой мощности, приходящейся на элемент, была определена ранее в выражениях (2.6), (2.9), (2.21) и (2.23). Мощность тепловыделения источников  из уравнения (2.4) задана в явном виде. Следовательно, все составляющие полной мощности, приходящейся на элемент на одном шаге по времени, вычислены. На следующем этапе расчета необходимо разработать метод интегрирования по времени.

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

Наиболее простым метод численного интегрирования обыкновенного дифференциального уравнения является метод Эйлера. Метод является явным, одношаговым, основан на представлении искомой функции в кусочно-линейном виде. Если представить уравнение теплового баланса в i-ого элемента (2.4) в виде

,

(2.27)

где

- полная мощность, приходящаяся на элемент на одном шаге по времени,

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

,

(2.28)

то метод Эйлер можно записать в виде

,

(2.29)

где

,  - температура i-ого элемента на n-ом и n+1-ом шагах по времени;

,  - теплоемкости i-ого элемента на n-ом и n+1-ом шагах по времени;

- задаваемый шаг по времени;

- полная мощность, приходящаяся на элемент на n-ом шаге по времени.

В случае учета зависимости теплоемкости от температуры, уравнение (2.29) решается методом Ньютона, в случае постоянной теплоемкости его можно переписать следующим образом

,

(2.30)

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

Очевидно, что обеспечить меньшую погрешность решения можно с использованием метода интегрирования более высокого порядка, чем метод Эйлера. Одним из наиболее широко используемых методов интегрирования является метод Рунге-Кутты четвертого порядка [], представляющий решение уравнения в виде полинома четвертой степени. Отличительной чертой метода является необходимость вычисления правой части уравнения (2.27) при трех дополнительных значениях температуры, распределенных между  и , что в случае проведения расчета радиационного теплообмена безматричным способом вынуждает проводить дополнительные циклы статистических испытаний. Это существенно увеличивает время расчета и является неприемлемым. Несколько менее распространенный метод Адамса [] позволяет получать решение с тем же порядком погрешности, что и метод Рунге-Кутты, но при этом не требует проводить дополнительные статистические расчеты. По методу Адамса решение строится на основе информации о нескольких предыдущих шагах интегрирования.

Следует отметить, что применение неявных методов интегрирования для решения задачи ограниченно теми же причинами, что и применение метода Рунге-Кутты. Неявные методы интегрирования связаны с решением системы нелинейных уравнений, правая часть которых не выражается аналитически и вычисляется методом Монте-Карло. При решении системы потребуется проводить дополнительные статистические расчеты правой части, что представляется излишне избыточным и ресурсозатратным.  

Таким образом, рациональным является интегрирование дифференциального уравнения (2.27) четырехшаговым явным методом Адамса с постоянным шагом, реализованным в следующем виде []:  

,

(2.31)

где

- полная мощность, приходящаяся на i-й элемент на n-ом, n−1-ом, n−2-ом, n−3-ом шагах по времени.

Аналогично методу Эйлера, в случае переменой теплоемкости уравнение (2.31) записывается в виде

,

(2.32)

и решается методом Ньютона. Когда теплоемкость принимается постоянной, температуру на следующем шаге можно получить в явном виде

.

(2.33)

На первом, втором и третьем шагах применить четырехшаговый метод Адамса не представляется возможным. Для расчета температуры на этих шагах необходимо применить метод Эйлера, двух- и трехшаговые методы Адамса соответственно []:

,

.

(2.34)

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

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

Рис. 2.14. Примеры решения задачи радиационно-кондуктивного теплообмена. Зависимости температур от времени

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

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

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

  •  рассчитываются входящие тепловые мощности для всех элементов модели,
  •  для каждого элемента вычисляются величины:

,

(2.35)

,

где

- шаг при охлаждении элемента,

- шаг при нагреве элемента,

,  - коэффициенты уменьшения и увеличения температуры за шаг соответственно,

- текущая температура на n−1-ом шаге по времени,

- полная мощность, приходящаяся на элемент на n−1-ом шаге по времени,

- удельная теплоемкость на n−1-ом шаге по времени,

- плотность,

- объем,

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

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

Рис. 2.15. Пример теоретической зависимости полной мощности элемента от времени

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

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

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

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

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

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

,

(2.36)

где

- скорректированная полная мощность для i-го элемента n-ом шаге,

- момент времени на начало n-ого шага,

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

,

,

,

,

(2.37)

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

,

(2.38)

где

- текущее значение шага интернирования, рассчитанное по формуле (2.35),

- скорректированное значение полной мощности для i-го элемента n-ом шаге интегрирования.

  1.   Генератор случайных чисел

Важнейшей составляющей метода Монте-Карло является механизм получения случайных чисел. Как было показано в [, ], результаты расчета по методу Монте-Карло прямым и непосредственным образом зависят от качества использованного генератора случайных чисел (ГСЧ), поэтому выбор подходящего ГСЧ является чрезвычайно важным этапом построения алгоритма.

Созданные к настоящему времени методы генерации истинно случайных чисел применяются в основном в криптографии для получения уникальных шифров. Непременным атрибутом ГСЧ является так называемый «источник энтропии», источник истинно случайного сигнала. В роли таких источников выступают счетчик тактов процессора, характеристики ЭВМ, действия пользователя, взаимодействия между запущенными на ЭВМ задачами и шумы звуковых карт. Большинству методов получения истинно случайных чисел присуща невысокая скорость работы и возможные сбои в генерации.

От многих недостатков ГСЧ избавлены программные генераторы псевдослучайных чисел (ГПСЧ). Их основой отличительной чертой является, то, что получаемые числа являются не совсем случайными, а имеют некоторый период, по истечении которого получаемые числа начинают повторяться.

Одним из первых ГПСЧ был неудачный алгоритм RANDU, результаты работы которого поставили под вопрос достоверность проведенных с его помощью исследований []. В настоящее время для различных целей применяется линейный конгруэнтный метод, берущий свое начало от алгоритма RANDU и широко используемый для получения псевдослучайных чисел в современных компиляторах. Также распространены линейный регистр сдвига с обратной связью и метод Фибоначчи с запаздываниями.

Наиболее широко применяется так называемый «вихрь Мерсенна», разработанный в 1996-1997 годах Макото Мацумото и Такудци Несимурой []. Метод основан на свойствах простых чисел Мерсенна и отличается высокой скоростью работы, равномерным распределением в 623 измерениях и периодом повторения, составляющим .

В настоящей работе для генерации псевдослучайных чисел используется вихрь Мерсенна в виде алгоритма MT19937.

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

Во второй главе настоящей работы описан метод расчета радиационно-кондуктивного теплообмена. Проверка достоверности и работоспособности метода было проведена путем оценки погрешности решения уже известных задач. Для этого метод был реализован в виде программы для ЭВМ.

  1.   Программная реализация метода

Метод, описанный во второй главе настоящей работы, реализован автором в программном модуле T.H.O.R.I.U.M. (Thermal optical radiation iteration module). Модуль предназначен для работы в среде операционной системы Microsoft Windows. T.H.O.R.I.U.M. написан на языке Visual Basic.NET. Код модуля является открытым и доступен в репозитории SourceForge.org по адресу http://www.sourceforge.net/projects/thorium.

Взаимодействие с OpenGL, обеспечивающим работу с аппаратными средствами графического ускорителя,  происходит посредством свободно распространяемой библиотеки классов Tao Framework 2.1.0.

Распределенные вычисления в локальной сети выполняются средствами операционной системы, в сети Интернет – с помощью протокола FTP.

Алгоритм вихря Мерсенна, используемого для генерации последовательности псевдослучайных чисел, так же реализован средствами языка Visual Basic.NET.

Среди особенностей реализации метода следует отметить следующее:

  •  совместимость по входным данным с широко распространенным комплексом конечно-элементного моделирования MSC.NASTRAN,
  •  представление результатов расчета как в текстовом виде, пригодном для последующего использовании в программных комплексах, так и в виде трехмерных графических изображений, пригодных для анализа человеком,
  •  возможность работы как в режиме с графическим интерфейсом, что позволяет пользователю проводить анализ входных данных на основе визуальной информации, так и в пакетном режиме с передачей аргументов расчета через командную строку Windows,  позволяющем встраивать T.H.O.R.I.U.M. в технологические цепочки проектирования изделий.

Входными данными модуля являются файл с конечно-элементной моделью конструкции, представленной в формате Bulk Data File комплекса MSC.NASTRAN [] и файл с описанием теплофизических и термооптических свойств материалов и покрытий, характеристиками источников излучения и тепловыделения, начальным полем температур и параметрами интегрирования.

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

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

Внешний вид T.H.O.R.I.U.M. в процессе расчета модели обсерватории «Миллиметрон» показан на Рис. 3.1.

Рис. 3.1. Внешний вид T.H.O.R.I.U.M. в процессе расчета модели обсерватории «Миллиметрон»

  1.   Подтверждение достоверности и работоспособности метода

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

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

3.2.1.   Радиационный теплообмен

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

Рис. 3.2. Две перпендикулярные прямоугольные пластины с общей стороной

Рис. 3.3. Две квадратные параллельные пластины

Зависимость угловых коэффициентов для перпендикулярных пластин выражается в следующем виде

,

(3.1)

где

, .


,

(3.2)

где

.

Для получения значений угловых коэффициентов в модуле T.H.O.R.I.U.M. для перпендикулярных пластин использовалось по 10000 пучков фотонов, излучаемых с каждой стороны пластины, для параллельных пластин  - по 1000 пучков с каждой стороны пластины. Аналитические и численные зависимости угловых коэффициентов перпендикулярных и параллельных пластин от их геометрических параметров изображены на Рис. 3.3-3.6.

Рис. 3.3. Аналитические и численные зависимости угловых коэффициентов перпендикулярных пластин .

Рис. 3.4. Аналитические и численные зависимости угловых коэффициентов перпендикулярных пластин .

Рис. 3.5. Аналитические и численные зависимости угловых коэффициентов перпендикулярных пластин

Рис. 3.6. Аналитические и численные зависимости угловых коэффициентов параллельных квадратных пластин

На Рис. 3.3-3.6 видно хорошее совпадение между аналитическими и численными результатами. Отсюда следует, что алгоритм расчета угловых коэффициентов реализован верно.

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

Аналогичная задача была решена в модуле T.H.O.R.I.U.M. с использованием метода интегрирования Эйлера с постоянным и переменным шагом, со сглаживанием и без сглаживания, а также с использованием четырехшагового метода Адамса.

Начальная температура обоих пластин составляла 300К, время интегрирования 100 с.

Таблица 3.1

Теплофизические свойства элементов модели

Материал

Плотность материала,

Удельная теплоемкость,

Теплопроводность,  

Толщина, мм

Лавсан с Al

1450

1000

0,24

1

Таблица 3.2

Термооптические свойства элементов модели

Коэффициент зеркального отражения

Коэффициент диффузного отражения

Коэффициент поглощения

Коэффициент пропускания

Степень черноты

0,0

0,0

1,0

0,0

1,0

Полученные зависимости температур и мощностей излучения от времени приведены на Рис. 3.7 и Рис. 3.8.

Рис. 3.7. Аналитическая и численные зависимости температуры от времени. Охлаждение параллельных пластин

Рис. 3.8. Аналитическая и численные зависимости мощности излучения от времени. Охлаждение параллельных пластин

На Рис. 3.7 и Рис. 3.8 видно, что аналитическая и численные зависимости температур и мощностей излучения от времени практически совпадают. Это позволяет сделать вывод о том, что описанные методы интегрирования реализованы корректно.

На следующем этапе проверялась правильность работы T.H.O.R.I.U.M. при расчете системы пластины из предыдущей проверки. В этом расчетном случае одна пластина поддерживалась при температуре 300К, а вторая пластина с начальной температурой 200К нагревалась от первой. Время интегрирования составляло 3000 с. Для расчета угловых коэффициентов с каждой пластин излучалось по 1000 пучков фотонов.

Полученные численные зависимости сравнивались зависимостями, полученными с помощью пакета MathCad. Результаты сравнения показаны на Рис. 3.9 и 3.10.

Рис. 3.9. Теоретическая и расчетная зависимости температуры от времени. Одна пластина поддерживается при постоянной температуре

Рис. 3.10. Теоретическая и расчетная зависимости мощности излучения от времени. Одна пластина поддерживается при постоянной температуре

Из анализа Рис. 3.9 и Рис. 3.10 следует, что теоретическая и расчетная температуры второй пластины совпадают с достаточной для практики точностью. Это позволяет утверждать, что 1000 пучков фотонов, излучаемых с каждой стороны пластины достаточно для получения приемлемых результатов. В том случае если необходимо повысить точность расчета, необходимо увеличить количество пучков, излучаемых с каждой стороны пластины.

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

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

Задача решалась с помощью пакета MathCad и модуля T.H.O.R.I.U.M. Результаты расчета представлены на Рис. 3.11-3.14.

Рис. 3.11. Теоретическая и расчетная зависимости температуры от времени. Освещаемая пластина

Рис. 3.12. Теоретическая и расчетная зависимости мощности излучения от времени. Освещаемая пластина

Рис. 3.13. Теоретическая и расчетная зависимости температуры от времени. Затененная пластина

Рис. 3.14. Теоретическая и расчетная зависимости мощности излучения от времени. Затененная пластина

Полученные результаты позволяют утверждать, что источник излучения в модуле T.H.O.R.I.U.M. реализован правильно.

Простейшей системой с переменными свойствами являются рассмотренные ранее параллельные пластины, выполненные из материала, теплоемкость которого зависит от температуры. Для проверки правильности расчета систем с переменной теплоемкостью эта задача была сначала решена с помощью MathCad, а затем с помощью модуля T.H.O.R.I.U.M.

Теплоемкость материала пластин зависела от температуры следующим образом, см. Рис. 3.15.

Рис. 3.15. Зависимость теплоемкости материала пластин от температуры

Постоянное значение теплоемкости соответствует приведенному выше расчету, в котором одна из пластин поддерживалась при постоянной температуре 300К, а другая, с начальной температурой 200К, нагревалась от нее.

Интерполированная по 20 точкам зависимость теплоемкости от температуры приведена на Рис. 3.15 для оценки точности интерполяции. Как было показано в разделе 2.4.2, в разработанном модуле T.H.O.R.I.U.M. для интерполирования таблично заданных функций используется метод обратных расстояний. Точность интерполяции по этому методу зависит от количества табулированных точек. На Рис 3.15 видно, что исходная линейная и интерполированная зависимости совпадают с достаточной степенью точности.

Результаты решения задачи с помощью MathCad  и T.H.O.R.I.U.M. приведены на Рис. 3.16 и Рис. 3.17. Анализ теоретического и расчетного решений позволяет сделать вывод об их близости, что, в свою очередь, позволяет утверждать, что расчет систем с переменной теплоемкостью реализован корректно.


Рис. 3.16. Теоретическая и расчетная зависимости температуры от времени. Теплоемкость зависит от температуры

Рис. 3.17. Теоретическая и расчетная зависимости мощности излучения от времени. Теплоемкость зависит от температуры

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

Зависимости коэффициентов диффузного отражения и поглощения от температуры показаны на Рис. 3.18 и Рис. 3.19.


Рис. 3.18. Зависимость коэффициента диффузного отражения от температуры

Рис. 3.19. Зависимость коэффициента поглощения от температуры

Интегрирование проводилось на интервале времени длиной 2000 с. На каждом шаге по времени с каждой из сторон пластин излучалось по 1000 пучков фотонов. Результаты расчета пластин с термооптическими коэффициентами, зависящими от температуры, приведены на Рис. 3.20 и Рис. 3.21.


Рис. 3.20. Теоретическая и расчетная зависимости температуры от времени. Термооптические коэффициенты зависит от температуры

Рис. 3.21. Теоретическая и расчетная зависимости мощности излучения от времени. Термооптические коэффициенты зависит от температуры

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

Зависимости коэффициентов диффузного отражения и поглощения от температуры и длины волны падающего излучения приведены на Рис. 3.22 и Рис. 3.23.

Интегрирование проводилось на интервале времени длиной 1000 с. На каждом шаге по времени с каждой из сторон пластин излучалось по 1000 пучков фотонов. Результаты расчета пластин с термооптическими коэффициентами, зависящими от температуры, приведены на Рис. 3.24 и Рис. 3.25.


Рис. 3.22. Зависимость коэффициента диффузного отражения от температуры и длины волны падающего излучения

Рис. 3.23. Зависимость коэффициента поглощения от температуры и длины волны падающего излучения

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

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

Анализируя Рис. 3.20, 3.21, 3.24 и 3.25 можно сделать вывод о том, что расчет радиационного теплообмена безматричным способом с учетом зависимости термооптических коэффициентов от температуры и длины волны реализован верно. Кроме того, следует отметить, что учет зависимости термооптических коэффициентов от длины волны позволяет существенно повысить точность расчетов, см. Рис. 3.24.

3.2.2. Кондуктивный теплообмен

Достоверность реализации расчета кондуктивного теплообмена проверялась на примере стержня и пластины с заданными источниками тепловыделения. Для этого одна и та же задача решалась с помощью программного пакета MathCad  и модуля T.H.O.R.I.U.M.

Исследуемый стержень имел длину 1 м и площадь поперечного сечения 0.0001 м. Плотность материала составляла 8960 . Расчет проводился для двух расчетных случаев. В первом случае теплоемкость и теплопроводность материала были приняты постоянными и составляли 386  и 401  соответственно. Во втором случае теплоемкость и теплопроводность зависели от температуры, см. Рис. 3.26 и Рис. 3.27. В обоих случаях радиационный теплообмен не учитывался, степень черноты поверхности стержня принималась равной 0,0.

Начальная температура стержня составляла 300К, на левом конце стержня находился источник тепла мощностью 1 Вт, на правом конце – охладитель мощностью -1 Вт.

В T.H.O.R.I.U.M. стержень моделировался 100 элементов типа CHEXA8.


Рис. 3.26. Зависимость теплоемкости материала стержня от температуры

Рис. 3.27. Зависимость теплопроводности материала стержня от температуры

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


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

Рис. 3.29. Теоретическая и расчетная зависимости мощности теплообмена по длине стержня от времени

Оценка достоверности расчета теплообмена в пластине проводилась при следующих исходных данных. Размеры пластины 1 м х 1 м, толщина пластины 0,01 м. Плотность материала пластины 8960 . Начальная температура пластины составляла 300К. В одном углу пластины находился источник тепла мощностью 25 Вт, в противоположном по диагонали углу пластины находился охладитель мощностью -25 Вт.  

Расчет пластины проводился для трех расчетных случаев: радиационный теплообмен не учитывается, постоянные теплоемкость и теплопроводность 386  и 401  соответственно; радиационный теплообмен не учитывается, теплоемкость и теплопроводность зависят от температуры, см. Рис. 3.26 и Рис. 3.27; учитывается радиационный теплообмен, теплоемкость и теплопроводность зависят от температуры, см. Рис. 3.26 и Рис. 3.27. В том случае, когда учитывается излучение с поверхности пластины, степень черноты принимается равной 1,0, в остальных случаях степень черноты 0,0.

Модель пластины в модуле состояла из 100 элементов типа CQUAD4.

Результаты расчета изображены на Рис. 3.30-3.32.


Рис. 3.30. Распределения температуры по пластине. Постоянные теплоемкость и теплопроводность. Нет излучения

Рис. 3.31. Распределения температуры по пластине. Теплоемкость и теплопроводность зависят от температуры. Нет излучения

На Рис. 3.30-3.32. видно, что температуры, полученные в пакете MathCad  и модуле T.H.O.R.I.U.M. достаточно хорошо совпадают. Это позволяет утверждать, что расчет кондуктивного теплообмена в конструкции, представленной двумерными элементами,  реализован корректно.

Рис. 3.32. Теоретическая и расчетная зависимости температуры по пластине от времени. Теплоемкость и теплопроводность зависят от температуры. Излучение учитывается

3.2.3. Заключение о достоверности разработанного метода расчета

В разделах 3.2.1. и 3.2.2. настоящей работы проведено исследование достоверности реализации и работоспособности модуля T.H.O.R.I.U.M. Исследование проводилось путем сравнения результатов, полученных с помощью разработанного метода, с известными аналитическими и численными решениями.

Оценивалась точность:

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

Результаты, полученные с помощью известных методов расчета и с помощью модуля T.H.O.R.I.U.M., с достаточной степенью точности согласуются между собой. Это позволяет утверждать, что метод расчета, разработанный в настоящей работе, является работоспособным и достоверным.

Достоверность и работоспособность метода позволяет рекомендовать модуль T.H.O.R.I.U.M. для практического использования.

  1.   Расчеты обсерватории «Миллиметрон»

Проект космической обсерватории «Миллиметрон» [] предполагает создание телескопа миллиметрового, субмиллиметрового и инфракрасного диапазона (10 мкм – 2 см) с диаметром главного зеркала 12 метров, охлажденным до температуры 4К. Обсерваторию планируется запустить на орбиту в районе точки L2 системы «Солнце-Земля». Ожидается, что «Миллиметрон» будет работать как в автономном режиме, так и в качестве космического плеча интерферометра с базами «Земля-космос». Проект «Миллиметрон» включен в Федеральную космическую программу России, запуск космического аппарата планируется в 2016 году. Внешний вид одного из вариантов обсерватории приведен на Рис. 1.11.

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

3.3.1. Оценка влияния погрешностей формы экранов на температурный режим зеркала

Пассивным средством обеспечения теплового режима обсерватории «Миллиметрон» являются радиационные экраны. Конструктивно экраны представляют собой относительно жесткий раскрывающийся каркас с закрепленной на нем тонкой металлизированной полиэтилентерефталатной пленкой типа лавсан (Mylar, Kapton). Толщина пленки составляет 10 мкм.

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

Для решения этой задачи были построены семь моделей системы «радиационный экран – зеркало» с различными отклонениями формы поверхности экрана.

Каждая модель состояла из трехслойного экрана и зеркала. Каждый слой экрана моделировался 100 элементов типа CQUAD4, зеркало моделировалось одним элементом того же типа. Расстояние между слоями экрана и между экраном и зеркалом составляло 0,1 м. Слои экрана были заданы параллельными друг другу и зеркалу. Отклонения формы задавались отдельно для каждого слоя экрана в 121 точке в направлении перпендикулярном плоскости экрана. Величина отклонения каждой из 121 точки поверхности выбиралась случайным образом из диапазона  , b назовем базой отклонения. Модели создавались с базами отклонения 0 мм, 1 мм, 2,5 мм, 5 мм, 7,5 мм, 10 мм, 20 мм. Моделируемая конструкция облучалась потоком солнечного излучения 225 со стороны экранов. Внешний вид модели с базой отклонения 20 мм показан Рис. 3.33.

Теплофизические и термооптические свойства экрана и зеркала приведены в Табл. 3.3 и Табл. 3.4 соответственно.

Рис. 3.33. Внешний вид модели с базой отклонения 20 мм

Табл. 3.3

Теплофизические свойства элементов модели

Элемент конструкции

Материал

Плотность материала,

Удельная тепло-емкость,

Тепло-проводность,  

Толщина, мм

Экран

лавсан

1450

1000

0,24

1

Зеркало

SiC

3200

750

120

1

Табл. 3.4

Термооптические свойства элементов модели

Поверхность

Коэф. зерк. отр.

Коэф. дифф. отр.

Коэф. погл.

Коэф. проп.

Степень черноты

Экран

0,0

0,95

0,05

0,0

0,05

Зеркало

0,0

0,95

0,05

0,0

0,05

Начальная температура всей системы  300К, время расчета  5.0E+5 с. Для интегрирования по времени использовался метод Эйлера с переменным шагом. Внешний точечный источник излучения мощностью 0.63E+26 Вт и эффективной температурой 5788К находился на расстоянии 1.5E+11 м от системы «зеркало-экраны».

Результаты расчета приведены на Рис. 3.34 и Рис. 3.35. Зависимости температуры зеркала от времени для различных значений базы отклонения b приведены на Рис. 3.34.


Рис. 3.34. Зависимость температуры зеркала от времени для различных значений базы отклонения b

Рис. 3.35. Зависимость температуры зеркала в момент времени 9.0E+5 от базы отклонения b

На Рис. 3.35 видно, что при увеличении базы отклонения температура зеркала увеличивается. При увеличении базы отклонения с 0 мм до 20 мм температуры зеркала возрастает с 55К до 95К.

Дополнительное исследование было проведено для проверки правильности приведенного расчета. Для этого были построены четыре модели.

Первая модель представляла собой один слой экрана с базой отклонения 0 мм без двух других слоев и без зеркала, см. Рис. 3.36. Вторая модель отличалась от первой базой отклонения, составлявшей 20 мм, см. Рис. 3.37. Все остальные параметры модели были приняты такими же, как и в предыдущем расчете в данном разделе, за исключением источника излучения. Его мощность составляла 3.85E+26 Вт, что обеспечивало лучистый поток на модель плотностью 1374 .

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


Рис. 3.36. Внешний вид модели одного слоя экрана с базой отклонения 0 мм

Рис. 3.37. Внешний вид модели одного слоя экрана с базой отклонения 20 мм

Высказанное предположение подтверждается расчетом: установившаяся средняя температура слоя с базой отклонения 0 мм составляет 331К, с базой отклонения 20 мм – 423К.

В третьей дополнительной модели внутренний, ближний к зеркалу слой экрана был выполнен с базой отклонения 20 мм, остальные слои – с базой с 0 мм.  В четвертой модели два внутренних слоя были выполнены с базой отклонения 20 мм. Плотность лучистого потока на модель составляла 225 . Рассчитывались установившиеся температуры зеркала и сравнивались с температурами зеркала в моделях со всеми слоями с отклонением 0 мм и 20 мм. Результаты сравнения представлены на Рис. 3.38.

Рис. 3.38. Зависимость температуры зеркала от числа искаженных слоев

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

3.3.2. Анализ влияния термооптических свойств покрытий на температуру зеркала

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

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

Табл. 3.5

Термооптические свойства элементов модели

Модель

Коэф. зерк. отр.

Коэф. дифф. отр.

Коэф. погл.

Коэф. проп.

Степень черноты

0,875

0,0

0,875

0,125

0,0

0,125

0,900

0,0

0,900

0,100

0,0

0,100

0,925

0,0

0,925

0,075

0,0

0,075

0,950

0,0

0,950

0,050

0,0

0,050

0,975

0,0

0,975

0,025

0,0

0,025

0,990

0,0

0,990

0,010

0,0

0,010

Каждая модель состояла из трехслойного экрана и зеркала сферической формы. Каждый слой экрана представлял собой шестигранную пирамиду с основанием, обращенным к зеркалу. Слой представлялся в модели 6 элементами типа CTRIA3. Зеркало моделировалось 92 элементами типа CQUAD4. Точечный источник излучения мощностью 3.85E+26 Вт и эффективной температурой 5788К находился на расстоянии 1.49E+11 м от зеркала на оси симметрии модели со стороны экранов. Источник обеспечивал лучистый поток плотностью 1374 на внешний слой экрана.

Размеры модели указаны на Рис. 3.39, внешний вид модели приведен на Рис. 3.40.

Рис. 3.39. Габаритные размеры модели

Рис. 3.40. Внешний вид модели в модуле T.H.O.R.I.U.M.

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


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

Рис. 3.42. Зависимость средней температуры зеркала в момент времени 1.0E+7 от коэффициента диффузного отражения

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

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

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

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

Модели экранов включали в свой состав источники излучения, обеспечивавшие лучистый поток плотностью 1374 на внешний экран. Во втором наборе моделей источники излучения отсутствовали.

Рассчитанные по этим моделям зависимости средних температур зеркала и внутреннего слоя экрана показаны на Рис. 3.43.

Рис. 3.43. Зависимости средних температур зеркала и внутреннего слоя экрана от коэффициента диффузного отражения

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

3.3.3. Сравнительный анализ двух вариантов конструкции обсерватории

Кроме варианта конструкции обсерватории, представленном на Рис. 1.11, Рис. 3.39 и Рис. 3.40, были предложены альтернативные варианты конструкции обсерватории. В частности, варианты конструкции с зеркалом диаметром 10 м, подвижно закрепленным относительно корпуса обсерватории.

Каждая модель состояла из шестислойного экрана и зеркала параболической формы с фокусным расстоянием 2,4 м. Каждый слой экрана представлял собой восьмигранную пирамиду с основанием, обращенным к зеркалу. Слой представлялся в модели восемью элементами. Зеркало моделировалось 108 элементами. Точечный источник излучения мощностью 3.85E+26 Вт и эффективной температурой 5788К находился на расстоянии 1.49E+11 м от зеркала на оси симметрии модели со стороны экранов и обеспечивал лучистый поток плотностью 1374 на внешний слой экрана. Начальная температура обсерватории составляла 300К, время расчета 1.0 E+7 с.

Внешний вид моделей описанных вариантов конструкции представлен на Рис. 3.44-3.47.


Рис. 3.44. Шаг экранов 1 м, зеркало на оси

Рис. 3.45. Шаг экранов 1 м, зеркало наклонено на 90

Рис. 3.46. Шаг экранов 0,2 м, зеркало на оси

Рис. 3.47. Шаг экранов 0,2 м, зеркало наклонено на 90

Представленные варианты отличаются друг от друга шагом между экранами, составляющем 1 м для варианта, изображенного на Рис. 3.44-3.45, и 0,2 метра для варианта на Рис. 3.49-3.47.

Габаритные размеры моделей показаны на Рис. 3.48 и Рис. 3.49.

Рис. 3.48. Габаритные размеры модели обсерватории

с шагом экранов 1 м. Размеры указаны в метрах

Рис. 3.49. Габаритные размеры модели обсерватории

с шагом экранов 0,2 м. Размеры указаны в метрах

Особенностью варианта с шагом экранов 0,2 метра являются двухслойные экраны. Каждый из шести слоев в этом варианте состоит из двух пленок, находящихся на расстоянии 10 мм друг от друга, см. Рис. 3.49.

Теплофизические свойства материалов обоих вариантов конструкции приведены в Табл. 3.6, термооптические свойства – в Табл. 3.7.

Табл. 3.6

Теплофизические свойства вариантов конструкции обсерватории

Элемент конструкции

Материал

Плотность материала,

Удельная тепло-емкость,

Тепло-проводность,  

Толщина, мм

Зеркало

Композит

С-С

1490

1130

0,45

5

Экраны

Лавсан

с Al

1450

1000

0,24

1

Табл. 3.7

Термооптические свойства вариантов конструкции обсерватории

Модель

Коэф. зерк. отр.

Коэф. дифф. отр.

Коэф. погл.

Коэф. проп.

Степень черноты

Рабочая сторона зеркала

0,0

0,900

0,100

0,0

0,100

Тыльная сторона зеркала

0,0

0,900

0,100

0,0

0,100

Освещаемый слой экрана

0,0

0,800

0,200

0,0

0,200

Следующий за освещаемым слой экрана

0,0

0,800

0,200

0,0

0,200

Прочие слои экрана

0,0

0,900

0,100

0,0

0,100

Результаты расчета температурных полей для различных вариантов приведены на Рис. 3.50 - 3.53. Зависимости средней температуры зеркала от времени для различных вариантов конструкции обсерватории приведены на Рис. 3.54.


Рис. 3.50. Шаг экранов 1 м, зеркало на оси. Температурное поле

Рис. 3.51. Шаг экранов 1 м, зеркало наклонено на 90. Температурное поле

Рис. 3.52. Шаг экранов 0,2 м, зеркало на оси. Температурное поле

Рис. 3.53. Шаг экранов 0,2 м, зеркало наклонено на 90. Температурное поле


Рис. 3.54. Зависимости средних температур зеркала от времени для различных вариантов конструкции обсерватории. Слева общий вид зависимостей, справа крупно показано время окончания расчета

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

3.3.4. Анализ влияния Земли на тепловой режим главного зеркала

При полете по орбите в районе точки Лагранжа №2 системы Земля-Солнце, обсерватория выходит из плоскости эклиптики на величину примерно 1,5 млн. км. Было высказано предположение о том, что существующая конфигурация экранов не обеспечивает одновременной защиты от излучения Солнца и от излучения Земли. Высказанное предположение подтвердилось. Выяснилось, что при неосевом освещении обсерватории Землей в течение 14 суток тепловой поток на зеркало увеличивается в два раза в сравнении с осевым освещением.

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

В расчете мощности излучения Земли принималось, что температура Земли составляет T=300К, степень черноты =0,7, площадь поверхности A=кв. м. Отсюда по закону Стефана-Больцмана получаем мощность излучения I:

Вт.

Расчетное положение Земли выбирается на основе анализа орбиты обсерватории, представленной на Рис. 3.55.

Рис. 3.55. Вид орбиты обсерватории во вращающейся системе координат.

Расстояния показаны в тыс. км

В инерциальной геоцентрической эклиптической системе координат ось Х в начальный момент времени направлена от Земли к Солнцу, ось Y лежит в плоскости эклиптики, ось Z – перпендикулярна плоскости эклиптики. В процессе полета продольная ось обсерватории поддерживается направленной на Солнце, при этом Земля совершает движение относительно этой оси.

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

.

Для проведения сравнения теплового режима был выбран и второй вариант расположения Земли. В этом случае ее координаты составляли:

.

Моделируемая обсерватория представляла собой экран, состоящий из 11 слоев и зеркала. Внешний вид модели обсерватории показан на Рис. 3.56, габаритные размеры модели показаны на Рис. 3.57.

Рис. 3.56. Внешний вид модели для исследования влияния Земли

Рис. 3.57. Габаритные размеры модели. Размеры указаны в метрах

Начальная температура зеркала обсерватории составляла 4,5К, ближайшего к зеркалу слоя экрана 42К, внешнего по отношению к зеркалу слоя – 320К. Температура остальных экранов линейно менялась от ближайшего к зеркалу до наиболее удаленного. Источниками излучения в системе выступали Солнце и Земля. В случае неосевого освещения Солнце находилось на расстоянии 150 млн. км. от внешнего экрана по оси телескопа, Земля на расстоянии 1 млн. км по оси телескопа от внешнего экрана и на расстоянии 1 млн. 430 тыс. км от оси обсерватории. В случае осевого освещения Солнце  находилось на том же расстоянии, Земля на расстоянии 1 млн. км по оси телескопа от внешнего экрана. Работа бортовых холодильных установок не учитывалась. Расчет проводился на 14 суток освещения. При моделировании были приняты следующие свойства теплофизические свойства материалов, Табл. 3.8.

Табл. 3.8

Теплофизические свойства вариантов конструкции обсерватории

Элемент конструкции

Плотность материала,

Удельная тепло-емкость,

Тепло-проводность,  

Толщина, мм

Зеркало

3200

750

120

5

Экраны

1450

1000

0,24

1

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

Расчетные зависимости максимальной температуры зеркала и мощности теплового потока на зеркало показаны на Рис. 3.58.

     

Рис. 3.58. Зависимости максимальной температуры зеркала и мощности теплового потока на зеркало для двух положений Земли

Анализируя зависимости, представленные на Рис. 3.58, можно сделать следующие выводы:

  •  максимальная температура зеркала по истечении 14 суток неосевого и осевого освещения составляет соответственно 4,8К и 4,653К.
  •  тепловой поток на зеркало по истечении 14 суток освещения в первом и втором случаях расположения Земли составляет соответственно 0,0686 Вт и 0,0340 Вт.

Видно, что проанализированная конфигурация экранов не обеспечивает теплового режима зеркала при заданном положении источников излучения.


Заключение

В заключении сформулированы основные результаты и выводы.

  1.  Предложен метод анализа радиационно-кондуктивного теплообмена в конструкции больших космических телескопов.
  2.  Предложен метод расчета лучистого теплообмена, опирающийся на подходы из области трехмерной графики.
  3.  Разработан алгоритм расчета больших космических телескопов, пригодный для оперативного анализа теплового режима на ЭВМ, обладающих ограниченными вычислительными ресурсами.
  4.  Разработан программный модуль T.H.O.R.I.U.M., реализующий предложенный алгоритм.
  5.  Получены результаты расчета тепловых режимов обсерватории «Миллиметрон»:
    1.  Показано существенное влияние погрешностей формы радиационного экрана на температуру зеркала телескопа.
    2.  Обнаружено существование набора термооптических коэффициентов, обеспечивающих наименьшую температуру зеркала при заданной геометрической конфигурации экранов.


Приложение. Краткое описание алгоритма T.H.O.R.I.U.M.

В настоящем документе приведено кратное описание алгоритма работы модуля T.H.O.R.I.U.M. с указанием основных вызываемых функций (Function) и процедур (Sub).

Приложение реализовано с помощью двух основных классов:

  •  Model, экземпляр, используемый в приложении, имеет имя FEM,
  •  TOP_Solver, экземпляр, используемый в приложении, имеет имя Solver.

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

П1. Подготовка расчета

Работа модуля начинается с запуска главной формы приложения Main_form (конструктор Main_form.New и событие Main_form.Load). В этих процедурах происходит начальная инициализация объектов FEM и Solver и других, вспомогательных внутренних переменных модуля, так же настройка вида главной формы. После этого на экране появляется форма, и приложение ждет действий пользователя. На этом этапе возможно:

  •  Выбрать файл с исходными данными, пункт меню File/Open;
  •  Выйти из приложения, пункт меню File/Exit;
  •  Получить справку о приложении, пункт меню About.

При выборе пункта меню File/Open запускается процедура считывания файлов исходных данных (Sub Main_form.Open_model).

П1.1. Считывание файлов с исходными данными

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

В случае выбора пользователем файла со скомпилированной моделью (файл с расширением *.model) процедура производит его десериализацию и заполняет поля объекта FEM.

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

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

П1.2. Отображение модели

Отображение модели производится по-разному в зависимости от степени завершенности расчета, записанного в файле.

Если пользователь выбрал файл с расширением *.bdf или *.dat, заведомо не содержащий результаты расчета, или файл с расширением *.model, не содержащий результаты расчета, то в этом случае создается экземпляр класса Visualisation, c именем view. Класс Visualisation используется для доступа к средствам OpenGL. Поля созданного объекта заполняются таким образом, чтобы на экране появилась вспомогательная форма с исходным видом модели. В данном режиме отображения возможен просмотр теплофизических свойств материалов элементов и термооптических свойств поверхностей.

Если пользователь выбрал файл с расширением *.model, в котором содержатся результаты расчета, то тогда на основной форме появляются объекты управления результатами, вспомогательная форма отображает модель, соответствующую выбранному шагу расчета.

После того как пользователь проверил свойства модели, он может запустить расчет выбором пункта меню Analysis/Thermal analysis.

П2. Проведение расчета

Переход к расчету производится после выбора пункта меню Analysis/Thermal analysis. С выбором этого меню скрывается вспомогательная форма для отображения модели view, так как она в дальнейшем в ходе расчета будет использоваться для определения взаимной видимости элементов. После того, как форма view стала невидимой, происходит переход к процедуре Solver.Solve_problem из класса TOP_solver.

П2.1. Подготовка объекта Solver к решению 

Процедура Solver.Solve_problem производит выбор того или иного алгоритма решения в зависимости от информации, содержащейся в файле с расширением *.top. В настоящем документе рассматривается проведение расчета с интегрированием по методу Эйлера с переменным шагом и учетом переменности теплофизических и термооптических свойств модели. Данный вид расчета реализуется процедурой TOP_Solver. Euler_var_T_step_with_smoothing_MATRIXLESS.

П2.2. Расчет методом Эйлера с переменным шагом 

Перед началом интегрирования по времени процедура создает отображение конечно-элементной модели в набор многоугольников-граней и присваивает всем гранями модели, принимающих участие в лучистом теплообмене (массив FEM.Shining_face), цвета по номерам. Расчет Соответствие номера цвету производится с помощью процедуры Face.Detect_face_number.

Созданный массив граней FEM.Polygon используется для создания экземпляра класса Visualisation с именем aux_Visual, который применяется для расчета попадания пучка фотонов в элемент.

Непосредственно перед началом интегрирования применяются начальные условия. Температуры и теплопритоки элементов приводятся в соответствие с заданными.

К интегрированию по времени все готово. Шаг интегрирования начинается с вызовов функций расчета кондуктивного теплообмена TOP_Solver. Calc_C_transfer и лучистого теплообмена TOP_Solver.Calc_R_transfer_N_photon.

П2.3. Расчет кондуктивного теплообмена 

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

П2.4. Расчет лучистого теплообмена 

В процедуре TOP_Solver.Calc_R_transfer_N_photon, реализующей расчет лучистого теплообмена методом Монте-Карло,  перед началом расчета обнуляются мощности поглощенного, зеркально и диффузно отраженного и пропущенного излучения (процедура Reset_fases_R_Q).

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

Имея рассчитанные мощности, становится возможным перейти к расчету мощности излучения, поглощаемого элементами. Для этого функцией TOP_Solver.Illuminate_model_and_save_absorbed_power_in_faces производится «освещение» модели с использованием в качестве источника излучения заданного элемента.

П2.5. «Освещение» модели

«Освещение» модели – это  расчет мощности излучения, испущенного заданным элементом и поглощенного всеми элементами модели. Эта задача решается с помощью функции TOP_Solver.Illuminate_model_and_save_absorbed_power_in_faces.

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

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

Аналогичным образом модели «освещается» из внешних источников излучения. Суммарная мощность теплопритока на i-й элемент на n-ом шаге по времени сохраняется в поле HT_element(i).HT_Step(n).Q, температура элемента на этом же шаге по времени хранится в поле HT_element(i).HT_Step(n).T.

П2.6. Сглаживание, выбор текущего шага и интегрирование мощности теплопритока

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

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

,

(1)

где

- скорректированная полная мощность для i-го элемента n-ом шаге,

- момент времени на начало n-ого шага,

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

,

,

,

,

(2)

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

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

  •  рассчитываются входящие тепловые мощности для всех элементов модели,
  •  для каждого элемента вычисляются величины:

,

(4)

,

где

- шаг при охлаждении элемента,

- шаг при нагреве элемента,

,  - коэффициенты уменьшения и увеличения температуры за шаг соответственно,

- текущая температура на n−1-ом шаге по времени,

- полная мощность, приходящаяся на элемент на n−1-ом шаге по времени,

- удельная теплоемкость на n−1-ом шаге по времени,

- плотность,

- объем.

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

Метод Эйлер с переменным шагом можно записать в виде

,

(5)

где

,  - температура i-ого элемента на n-ом и n+1-ом шагах по времени;

,  - теплоемкости i-ого элемента на n-ом и n+1-ом шагах по времени;

- задаваемый шаг по времени;

- полная мощность, приходящаяся на элемент на n-ом шаге по времени.

В случае учета зависимости теплоемкости от температуры, уравнение (5) решается методом Ньютона, в случае постоянной теплоемкости его можно переписать следующим образом

,

(6)

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

Уравнение (5) решается до достижения заданного времени окончания интегрирования. После этого приложение переключается в режим отображения результатов.

П3. Отображение и сохранение результатов расчета

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

П3.1.Отображение результатов расчета

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

П3.2.Сохранение результатов расчета

Сохранение результатов возможно в двух видах: в виде сериализованной модели в файле с расширением *.model и в виде таблицы температуры и мощности теплопритока элементов в текстовом файле *.txt. Первый вид сохранения вызывается пользователем при выборе пункта меню File/Save, второй при выборе пункта меню File/Save results as txt. В первом случае вызывается процедура Main_form.Save_model, с помощью диалогового окна запрашивающая у пользователя имя файла и сериализующая объекты FEM и Solver в указанный файл. Во втором случае вызывается процедура Main_form.Save_txt_results, сохраняющая в выбираемый пользователем файл массив HT_element(i).HT_Step в виде набора строк: «Номер_элемента Время Температура Мощность_теплопритока».

П4. Распределенный параллельный расчет

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

Подготовительный этап

  1.  Выделяется ЭВМ-сервер и задается ее адрес в локальной сети или Интернет.
  2.  На сервер устанавливается программное обеспечение в серверном варианте для проведения распределенных вычислений.
  3.  На ЭВМ-клиенты устанавливается и запускается программное обеспечение в клиентском варианте, значение поля Main_form.DC_Mode = "CLIENT".

Этап проведения расчета

  1.  На сервер по заранее заданному адресу помещается модель конструкции, тепловой расчет которой необходимо провести.
  2.  Происходит запуск приложения в серверном варианте, значение поля Main_form.DC_Mode = "SERVER".
  3.  Сервер производит считывание модели, процедура Main_form.Open_model,  и сериализует ее по указанному адресу Main_form.Save_model.
  4.  Сервер формирует вектор температур на текущем шаге времени с помощью процедуры TOP_Solver.Write_temperature_file_to_dir.
  5.  Клиенты через заданный промежуток времени проверяют наличие модели на сервере и при наличии модели производят ее считывание, процедура TOP_Solver.Read_temperature_file_from_dir.
  6.  Клиенты производят оценку собственной производительности функцией TOP_Solver. Calc_productivity_index и сообщают о ее результатах серверу.
  7.  Сервер, получив результаты оценки производительности, составляет списки элементов, процедура TOP_Solver.Detect_part_of_model_for_clients, лучистый поток с которых необходимо рассчитать клиентам, и передает эти списки клиентам. Количество элементов, которые необходимо рассчитать клиенту, прямо пропорционально его производительности: чем больше вычислительная мощность клиента, тем больше элементов ему необходимо рассчитать.
  8.  Сервер после рассылки списков для расчета лучистого теплообмена самостоятельно проводит расчет кондуктивного теплообмена TOP_Solver.Calc_C_transfer.
  9.  Клиенты, производят расчет мощности излучения, испускаемого элементами из списка и поглощаемого всеми элементами модели TOP_Solver.Calc_R_transfer_N_photon. По окончании расчета клиент формирует вектор поглощенных мощностей и передает его на сервер TOP_Solver.Write_heat_power_file_to_dir.
  10.  Клиенты после отправки вектора поглощенных мощностей на сервер переходят к пункту 5.
  11.  Сервер, закончив проведение расчета кондуктивного теплообмена, ожидает окончания расчетов клиентов. В случае, если по истечении заданного срока все клиенты не предоставили векторы поглощенных мощностей излучения, то сервер производит недостающие расчеты самостоятельно.
  12.  Сервер, получив все векторы поглощенных мощностей излучения, суммирует их, определяет суммарные тепловые потоки от лучистого и кондуктивного теплообмена и производит расчет вектора температур на следующем шаге времени TOP_Solver.Euler_var_MATRIXLESS_DC.
  13.  Сервер размещает новый вектор температур и переходит к пункту 6.
  14.  По окончании расчета сервер формирует по заданному адресу сигнал клиентам об окончании расчета. При получении этого сигнала клиенты переходят в режим ожидания расчета, увеличивая интервал времени между проверками наличия модели на сервере.

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


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

  1.  

ABAQUS Analysis User’s Manual v.6.5. – Velizy-Vallacoublay: Dassault Systemes S.A., 2004. – 1467 p.

  1.  

Anderson G. E. Thermal Radiation Analyzer System (Cray Version with NASADIG) // NASA Technical Report. – 1994, Document ID: 19940002537.

  1.  

ANSYS Thermal Analysis Guide v.10.0. – Canonburg: Ansys Inc., 2005. – 82 p.

  1.  

Arvo J. et al. Monte Carlo Ray Tracing // Proc. Siggraph 2003. – 2003. Cource 44. 171 p.

  1.  

Ashdown I. Eigenvector Radiosity, PhD thesis. –Vancouver: Department of Computer Science, The University of British Columbia, 2001. – 138 p.

  1.  

Bamberg J.A., Zaun N.H. Design and performance of the cryogenic focal plane optics assembly for the Infrared Astronomical Satellite (IRAS) // Proc. Cryogenic optical systems and instruments. – 1984. – P. 94-102

  1.  

Baryshev A.M., Kardashov N.S., Arkhipov M.Yu. et al. Deployable Antennas for Space Radio Telescope: Radioastron and Millimetron Missions // Proceedings of the 30th ESA Antenna Workshop. – 2008.

  1.  

Behee R. Introduction to Thermal Modeling Using SINDA/G: A Tutorial Guide. –Chandler: Network Analysis Inc., 2006. –157 p.

  1.  

Bennett C.L. et al. The Microwave Anisotropy Probe (MAP) Mission // Astrophysical Journal. – 2003. – Vol. 143 P. 567-576

  1.  

Brewster M. Q. Thermal Radiative Transfer and Properties. – New York: John Wiley & Sons, Inc., 1992. – 356 p.

  1.  

Burrows C.J., Holtzman J.A., Faber S.M., Bely P.Y., Hasan H., Lynds C.R., Schroeder D. The imaging performance of the Hubble Space Telescope // Astrophysical Journal. – 1991.—Vol.369. –P. L21-L25

  1.  

Catmull E.E. A subdivision algorithm for computer display of curved surfaces. – Salt Lake City: Univercity of Utah, 1974. – 83 p.

  1.  

Cohen M.F., Chen S.E., Wallace J.R., Greenberg D. P. A progressive refinement approach to fast radiosity image generation // Proc. of SIGGRAPH. – 1988. – Vol. 22.–P.75-84

  1.  

Cohen M.F., Greenberg D.P. Hemi-Cube: A Radiosity Solution For Complex Environments // Proc. Of the 12-th annual conference on Computer Graphics and interactive techniques. – 1985. –Vol. 19. –P.31-40

  1.  

Cook R.L., Porter T., Carpenter L. Distributed ray tracing // Proc. of SIGGRAPH. – 1984. – Vol. 18.–P.137-145

  1.  

DiPirro M. et al. Cooling Technology for Large Space Telescopes // Proceedings of the SPIE. – 2007. Vol. 6687. – P. 66870D-1-66870D-11

  1.  

DoE Fundamentals handbook. Thermodynamics, Heat Transfer, and Fluid Flow, Vol 2., DOE-HDBK-1012/2-92, 1992, 58 p.

  1.  

Eckart P. The Lunar Base Handbook. – New-York: McGraw-Hill, 2006. – 820 p.

  1.  

Efstathiou G., Lawrence C., Tauber J. Plank. The scientific programme. – ESA-SCI(2005)1, 2005. 152 p.

  1.  

Fazio G.G. Recent Results from the Spitzer Space telescope: A New View of the Infrared Universe // Neutrinos and Explosive Events in the Universe. – 2005. – Vol. 209. P. 47-71  

  1.  

Frisk U.O. Status of the ODIN Project // Infrared and Submillimeter Space Astronomy. – 2002. – Vol. 4 – P. 211-217

  1.  

Fuchs H., Kedem Z.M., Naylor B.F. On visible surface generation by a priori tree structures // Proc. of SIGGRAPH. – 1980. – Vol. 14.–P.124-133

  1.  

Gardner J.P. et al. The James Webb Space Telescope // Space Science Reviews. – 2006. – Vol. 123. – P. 485-606

  1.  

Gilmore D.G., Bello M. Satellite Thermal Control Handbook. – El Segundo: Aerospace Corporation, 2002. – 581 p.

  1.  

Glassner A.S. Space subdivision for fast ray tracing // IEEE Computer Graphics & Application. – 1984. – Vol. 4.–P.15-22

  1.  

Goble R.G., Jensen C.L. Thermal radiation analysis system TRASYS 2: User’s manual Revision 3. – Denver: Martin Marietta Corp., 1980. – 100 p.

  1.  

Golliher E. A comparison of TSS ans TRASYS in form factor calculation // Proc. 5th Annual Thermal and Fluids Analysis Workshop.—1993. – P.41-46

  1.  

Goodwin P.S., Meeks W. G., Morris R. E. Helios mission support // The Deep Space Network. – 1976. – P. 31-37

  1.  

Goral C., Torrance K.E., Greenberg D.P., Battalie B. Modeling the interaction of light between diffuse surfaces // Proc. of SIGGRAPH. – 1984. – Vol. 18.–P.213-222

  1.  

Gorenstein P., Harris B., Gursky H., Giacconi R. A rocket payload using focusing X-ray optics for the observation of soft cosmic X-rays // Nuclear Instruments and Methods. – 1971. – Vol. 91. – P. 451-459

  1.  

Gursky, H., Schnopper, H., Parsignault, D. The Hard X-ray experiment on the Astronomical Netherlands Satellite // Astrophysical Journal. – 1975 –  Vol. 201. – P. L127-L131

  1.  

Harris M. Mapping computational concepts to GPUs // ACM SIGGRAPH 2005 Cources. – 2005.

  1.  

Harwit M. The Herschel Mission // Advances if Space Research. – 2004. – Vol. 34. – P. 568-572

  1.  

Heller C. Interface Software Development for Patran/Thermal, Esarad and Thermica // Proc. 19th European Workshop on Thermal and ECLS Software.—2005. – P.K-1-K-15

  1.  

Ishihara D. et al. Mid-infrared all-sky survey with AKARI // Memorie della Societa Astronomica Italiana. – 2006. – Vol. 77. – P. 1089-1094

  1.  

Jackson C. Integration of Monte Carlo Methods with TMG’s Suite of Radiation Analysis Tools // Proc. 19th European Workshop on Thermal and ECLS Software.—2005. – P.H-1-H-12

  1.  

Jensen H.J., Christensen N.J. Photon maps in bidirectional monte carlo ray tracing of complex objects // Computer Graphics. – 1995. – Vol. 19.–P.215-224

  1.  

Kessler M.F. et al. The ISO mission – A scientific Overview. ESA Bulletin № 84, 1995. – 43 p.

  1.  

Kirtley C. ESATAN Thermal suite status // Proc. 20th European Workshop on Thermal and ECLS Software.—2006. – P.133-146

  1.  

Lawson D. TPF-C Technology Plan. Pasadena: Jet Propulsion Laboratory, 2005. – 154 p.

  1.  

Lienhard IV J.H., Lienhard V J.H. A Heat Transfer Textbook. – Cambridge, Phlogiston Press, 2003. 750 p.

  1.  

Liszka T. An interpolation method for an irregular net of nodes // International Journal for Numerical Methods in Engineering. – 1983.—Vol. 20. – P.1599 -1612

  1.  

Long C., Sayma N. Heat Transfer. – Holstebro: Ventus, – 2006. – 156 p.

  1.  

Love T.J. Radiative Heat Transfer. – Charles E. Merryl Books, Inc., 1968. – 118 p.

  1.  

Lunar Reconnaissance Orbiter Overview: The instrument Suite and Mission // Space Science Reviews. – 2007. – Vol. 129. – P. 391-419

  1.  

M. Reymond. MSC.Nastran 2005 DMAP Programmer’s Guide, Santa Ana, MSC.Software, 2005, Vol. 1. 1734 p.

  1.  

Mahan J. R.  Radiation Heat Transfer: A Statistical Approach, New York, Wiley, 2002, 504 p.

  1.  

Malley T. A shading method for computer generated images. Master’s thesis. – Salt Lake City: Univercity of Utah, 1988. – 111 p.

  1.  

Marsaglia G. Random numbers Fall Mainly in the Planes // Proc. of National Academy of Sciences. – 1968. – Vol. 61. –P.25-28

  1.  

Matsumoto М., Nishimura T., Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator // ACM Transactions on Modeling and Computer Simulation. – 1998. – Vol. 8. , – P. 3-30

  1.  

Metropolis N., Ulam S. The Monte Carlo Method // J. Amer. Statistical association. – 1949. Vol.  44. – P. 335-341

  1.  

Modest M.F. Radiative Heat Transfer. – San Diego: Academic Press, 2003. –

822 p.

  1.  

Muzumder S., Kersch A. A fast Monte Carlo Scheme for Thermal Radiation in Semiconductor Processing Application // Numerical Heat Transfer, Part B: Fundamentals. – 2000. – P. 185-199

  1.  

NEVADA (Net Energy Verification and Determination Analyzer). RENO and VEGAS. Quick Reference Guide – Incline Village: TAC Technologies, 2000. – 25 p.

  1.  

Nusselt W. Graphische bestimmung des winkelverhaltnisses bei der wärmestrahlung // Zeitschrift des Vereines Deutscher Ingenieure. – 1928. –Vol. 70. – P.673

  1.  

Odenwald S.J., Newmark J., Smoot G. A study of external galaxies detected by COBE Diffuse Infrared Background Experiment // Astrophysical Journal. – 1998. – Vol. 500. – P. 554-568

  1.  

Owens J.D., Luebke D. et al. A Survey of General-Purpose Computation on Graphics Hardware // Computer Graphics Forum. – 2007. –Vol. 26. – P.80-133

  1.  

Paige D.A. et al. The Lunar Reconnaissance Orbiter Diviner Lunar Radiometer Experiment // Space Science Reviews. – 2010. – Vol. 150. – P. 125-160

  1.  

Panczak T. Non-Grey and Temperature Dependent Radiation Analysis Methods // Proc. Thermal & Fluids Analysis Workshop 2005. – 2005. Short Cource

  1.  

Revnivtsev M.G., Sunyaev R. A., Gilfanov M. R., Churazov E.M.,

Goldwur A., Paul J., Mandrou P., Roques J. P. A Hard X-ray Sky Survey with the SIGMA Telescope of the GRANAT Observatory // Astronomy Letters. – 2004. –Vol. 30. P. 527–533

  1.  

Richer S.W. Experimental Determination of In Situ Utilization of Lunar Regolith for Thermal Energy Storage // NASA Contractor Report 191050. – 1993. – P. 1-10

  1.  

Rubinstein R.Y., Kroesse D.P. Simulation and the Monte Carlo Method. – New York: John Wiley & Sons, Inc., 2007. – 372 p.

  1.  

Schatz M. et al. High-throughput sequence alignment using Graphics Processing Units // BMC Bioinformatics. – 2007. – Vol.8. – P.1-10

  1.  

Shaughnessy B.M., Newborough M. A New Method for Tracking Radiative Paths in Monte Carlo Simulation // ASME Journal of Heat Transfer. – 1998. – Vol. 120. – P. 792-795

  1.  

Shaughnessy B.M., Newborough M. Calculating Reflected Paths of Radiation in High Reflectivity Enclosures // Proc. 5th ASME/JSME Joint Thermal Engineering Conference. – 1999

  1.  

Shepard D. A two-dimensional interpolation function for irregularyly-space data // Proc. Of the 1968 23rd ACM national conference.—1968.—P.517-524 

  1.  

Shirley P., Wang C., Zimmerman K.. Monte carlo techniques for direct

lighting calculations // ACM Trans. Graph.. – 1996. – Vol. 15. P. 1–36

  1.  

Shreiner D. OpenGL reference manual: the official reference document to OpenGL, version 1.4 – Boston: Addison-Wesley, 2004. – 760 p.

  1.  

Soriano T. Review of Improvements in Thermica // Proc. 19th European Workshop on Thermal and ECLS Software.—2005. – P.K-1-K-15

  1.  

Sparrow E.M., Cess R.D. Radiation Heat Transfer. – Brooks/Cole, 1978. – 185 p.

  1.  

Thermal Synthesizer System User manual. Volume 3. v.12.01. –  Houston: Spacedesign Corporation, 2007. – 55 p.

  1.  

THERMICA User’s Manual v.3.2. – Toulouse: Astrium, 2003. – 504 p.

  1.  

Thomas J. ESARAD status // Proc. 20th European Workshop on Thermal and ECLS Software.—2006. – P.79-92

  1.  

Vogt R.A. TRASYS – Thermal Radiation Analyser System (DEC VAX Version without NASADIG) // NASA Technical Report Server 19940003065. – 1994. – P. 1-10

  1.  

Vueghs P.E. Innovative Ray Tracing Algorithms for Space Thermal Analysis. – PhD thesis, Department of Aerospace and Mechanics, Liege, The University of Liege, 2009, 234 p.

  1.  

Wallace J.R., Elmquist K.A., Haines E.A. A Ray Tracing Algorithm for Progressive Radiosity // Proc. of SIGGRAPH. – 1989. – Vol. 23.–P.315-324

  1.  

Whitted T. An improved illumination model for shaded display structures // Communication of the ACM. – 1980. – Vol. 23.–P.343-349

  1.  

Wiebelt J. A. Engineering Radiation Heat Transfer. – Holt, Rinehart and Winston, New York, 1966. – 278 p.

  1.  

Wild W., Kardashev N. S. et al. Millimetron–a large Russian-European submillimeter space observatory // Experimental Astronomy. – 2008. – Vol. 23. – P. 221-244

  1.  

Wright R.S., Lipchak B. OpenGL SuperBible. – Indianapolis: Sams, 2004. – 1173 p.

  1.  

Ziering M. B., Sarofim A. F.. The electrical network analog to radiative

transfer: Allowance for specular reflection // ASME Journal of Heat Transfer. – 1966. – Vol. 88. – P. 341–342

  1.  

Авдуевский В.Д., и др. Основы теплопередачи в авиационной и ракетно-космической технике. М.: Машиностроение, 1992. – 528 с.

  1.  

Адинец А.В. Высокоуровневая система программирования графических процессорных устройств. Автореферат на соискание ученой степени кандидата физико-математических наук. – М.: МГУ им. М.В. Ломоносова, 2009. – 19 с.

  1.  

Алексеев В. П., Вайнштейн Г.Е., Герасимов П.В. Расчет и моделирование аппаратов криогенных установок. Л.: Энергоатомиздат, 1987. – 278 с.

  1.  

Блох А.Г., Журавлев Ю.А. Рыжков Л.Н. Теплообмен излучением. –М.: Энергоатомиздат, 1991. – 432 с.

  1.  

Виноградов И.С. Радиационное охлаждение зеркала крупногабаритного

космического телескопа. В кн. Радиоастрономическая техника и методы. М.: ФИАН, 2000 – Т.228, C. 112-128

  1.  

Виноградов И.С., Архипов М.Ю. Моделирование теплового и термодеформированного состояния рефлектора криотелескопа.

Этап 1: Выбор и уточнение методики, разработка модели, подготовка исходных данных для расчета температурных полей // Труды Российской конференции пользователей MSC 2006. – 2006 – с. 1-9

  1.  

Ву М., Девис Т., Нейдер Дж., Шрайнер Д. OpenGL. Руководство по программированию. Библиотека программиста.—СПб: Питер, 2006. – 623 с.

  1.  

Дульнев Г.Н., Парфенов В.Г., Сигалов А.В. Применение ЭВМ для решения задач теплообмена. – М.: Высшая школа, 1990. – 207 с.

  1.  

Елисеев В.Н., Белоногов Е.К., Расчеты тепловых режимов с использованием ЭВМ. Учебное пособие по курсу «Тепловые режимы космических объектов». – М.: МВТУ им. Н.Э. Баумана, 1989. – 47 с.

  1.  

Залетаев В.М., Капинос Ю.В., Сургучев О.В. Расчет теплообмена космического аппарата. – М.: Машиностроение, 1971. –  208 с.

  1.  

Землянский Б.А. Пакет программ ТЕРМ. Описание программы. – Королев: ЦНИИМаш, 2006. – 105 с.

  1.  

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

  1.  

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

  1.  

Кнут Д. Искусство программирования. Т.2 Получисленные алгоритмы. – М.: Вильямс, 2006. – 832 с.

  1.  

Крейт Ф., Блэк У. Основы теплопередачи. – М.: Мир, 1983. – 512 с.

  1.  

Лыков А.В. Тепломассообмен. Справочник. – М.: Энергия, 1972. – 560 с.

  1.  

Оцисик М.Н. Сложный теплообмен. М.: Мир, 1976. – 615 с.

  1.  

Планк М., Теория теплового излучения. – М.: Гостехтеориздат, 1935. – 209 с.

  1.  

Самарский А.А., Гулин А.В. Численные методы. – М.: Наука, 1989. – 432 с.

  1.  

Херн Д., Бейкер М.П. Компьютерная графика и стандарт OpenGL.—М.: Вильямс, 2005. – 1158 с.

  1.  

Цветков Ф.Ф., Григорьев Б.А., Тепломассообмен: Учебное пособие для вузов. – М.: МЭИ, 2005. – 550 с.

  1.  

Шикин Е.В., Боресков А. В. «Компьютерная графика. Полигональные модели», Москва, Диалог-МИФИ, 2000. – 464 с.


 

А также другие работы, которые могут Вас заинтересовать

62356. ТРЕБОВАНИЯ К УРОКУ РУССКОГО ЯЗЫКА В НАЧАЛЬНЫХ КЛАССАХ 33.42 KB
  Целенаправленность чёткость и логика урока Вопросы для анализа самоанализа урока: Какие цели поставлены учителем на уроке Соответствуют ли они теме и требованиям учебной программы Цели сформулированные учителем при подготовке к уроку должны соответствовать теме урока и требованиям учебной программы.
62357. Тренировочные упражнения к уроку русского языка в 11 классе 20.67 KB
  Много ошибок допускается при работе с предложением такого типа: Опровергая общее утверждение 1 приведите хотя бы один довод против. А в варианте 1 ответа предложение определённо-личное то есть подлежащего нет но сказуемое глагол 2ого лица множественного числа косвенно его выражает опровергая приведите довод против.
62358. Робота з природним матеріалом. Флоромозаїка з насіння 1.93 MB
  Мета: дати учням поняття про флоромозаїку; розширити знання про історію виникнення мозаїки; навчити дітей виготовляти композицію з насіння: Квітка осені сприяти розвитку уяви фантазії мислення уваги естетичного смаку образного мислення; мілкої моторики мязів пальців; сприяти вихованню охайності уважності та інтересу до роботи.
62361. Вправи та задачі на засвоєння таблиць додавання і віднімання. Кругові вирази. Побудова геометричних фігур за зразком 594.61 KB
  Діти до нас сьогодні в гості завітав космонавт. Лист Любі діти я тільки прилетів з дуже далекого космосу. На яке число ми щойно збільшували і зменшували інші числа на 6 Діти відкрийте свої зошити знайдіть новий робочий рядок і запишіть каліграфічно це число в рядок через клітинку.
62364. Слова, які означають назви предметів (іменники) 354.96 KB
  До кожного з них поставте запитання хто або що вчитель читає слова а учні хором ставлять запитання хто або що Комбайн комбайнер льотчик літак футбол село селянин корівник корова мурашка мурашник. На яке питання відповідає слово відгадка Що означає це слово назву предмета...