Содержание к диссертации
Введение
ГЛАВА 1. Обзор работ по методу дискретных ординат 18
ГЛАВА 2. Метод дискретных ординат на ортогональных сетках 30
2.1. Уравнение переноса излучения 30
Трехмерная геометрия 32
Двумерная геометрия 33
Цилиндрическая геометрия 33
Сферическая геометрия 35
2.2 Метод дискретных ординат на ортогональных сетках, основные расчетные алгоритмы 36
Трехмерная геометрия 37
Двумерная геометрия 43
Цилиндрическая геометрия 45
ГЛАВА 3. Выбор весовых функций 54
3.1 Плоская геометрия 54
3.2 Двумерный случай 58
3.3 Трехмерная геометрия 64
3.4 Цилиндрическая геометрия 64
ГЛАВА 4. Метод дискретных ординат на треугольных (тетраэдальных) неструктурированных сетках 66
Двумерная геометрия 66
Трехмерная геометрия 70
Алгоритм расчета 74
Расчетные сетки 78
ГЛАВА 5. Метод дискретных направлений 81
ГЛАВА 6. Квадратуры на поверхности сферы 85
ГЛАВА 7. Применение метода дискретных ординат на различных типах сеток 94
Двумерная геометрия 94
Цилиндрическая геометрия 95
Результаты тестирования разработанных программных кодов 101
ГЛАВА 8. Расчет радиационных потоков на внутренюю поверхность водородного и воздушного лазерного плазменного генератора 104
8.1 Постановка задачи радиационной газовой динамики ЛПГ 105
8.2 Водородный лазерный плазменный генератор 108
Свойства водородной плазмы 108
Результаты численного моделирования водородного ЛПГ . 109
8.3 Воздушный лазерный плазменный генератор 121
Свойства воздушной плазмы 121
Результаты численного моделирования воздушного ЛПГ 124
ГЛАВА 9. Расчет радиационных потоков на поверхность спускаемого космического аппарата при входе в атмосферу марса 129
Результаты вычислений 129
Заключение 137
Литература 138
- Метод дискретных ординат на ортогональных сетках, основные расчетные алгоритмы
- Цилиндрическая геометрия
- Результаты тестирования разработанных программных кодов
- Результаты численного моделирования водородного ЛПГ
Введение к работе
Теплообмен излучением играет важную роль в природе и технике. Структура атмосфер планет и звездных фотосфер, рабочие процессы в камерах сгорания и плазменных генераторов, тепловой режим радиоэлектронной аппаратуры и искусственных спутников Земли, аэротермодинамика космических аппаратов- это лишь некоторые примеры процессов, в которых теплообмен излучением является определяющим. В последнее время в научной литературе по сложному теплообмену отмечается повышенный интерес к радиационному переносу энергии в связи с его принципиальным значением для таких объектов техники, как возвращаемые космические аппараты, различные энергетические установки и термоядерные устройства.
Использование космических аппаратов для исследования планет солнечной системы, вызывает практический и теоретический интерес к изучению процессов протекающих при вхождении спускаемых аппаратов в плотные слои атмосфер планет. При скоростях составляющих 6-11 километров в секунду космический аппарат подвергается сильному конвективному и радиационному нагреву в результате взаимодействия с газом, нагретым в ударном слое.
Чтобы уменьшить общий вес космического аппарата необходимо точно рассчитать его тепловую защиту. Для этого важно с высокой достоверностью вычислить тепловые потоки на поверхность спускаемого модуля. Как правило, оценку конвективного и радиационного нагрева проводят для головного аэродинамического щита, принимающего на себя основную часть тепловой нагрузки. Однако при входе космического аппарата в атмосферу Марса радиационный нагрев задней поверхности может играть важную роль, поскольку эта поверхность практически не защищена, но подвергается воздействию радиационного теплового потока, испускаемого десятками кубических метров нагретой до высоких температур двуокиси углерода, которая, как известно, является хорошим излучателем в инфракрасной области спектра. Таким образом, одной из задач аэротермодинамического анализа межпланетных спускаемых аппаратов является предсказание интенсивности радиационного нагрева задней поверхности.
Представляет значительный практический и научный интерес создание теоретических моделей, разработка численных методов и программ расчета спектральных и интегральных радиационных тепловых потоков к поверхности космических аппаратов сложной формы и энергетических устройств типа плазменных генераторов. Это выдвигает необходимость разработки двух- и трехмерных моделей радиационного переноса тепла в объемах сложных геометрий на структурированных и неструктурированных сетках. Цель работы и задачи исследования.
Диссертационная работа имеет четыре основных цели:
Разработка численных методов для расчета спектральных и интегральных параметров поля излучения низкотемпературной плазмы в объемах сложной геометрии.
Разработка численных методов расчета излучательной способности локализованных в пространстве неоднородных объемов низкотемпературной плазмы.
Для этих целей были решены следующие задачи:
Разработаны алгоритмы решения уравнения переноса излучения на ортогональных сетках методом дискретных ординат в многогрупповом приближении, в плоской и цилиндрической геометриях;
Разработаны алгоритмы решения уравнения переноса излучения на неструктурированных треугольных и тетраэдальных сетках методом дискретных ординат в многогрупповом приближении;
Разработан алгоритм оптимального счета на неструктурированных сетках, а так же алгоритм расчета геометрических свойств ячеек и распознавания их пространственной ориентации относительно направления распространения фотонов;
Рассчитаны угловые квадратуры высоких порядков для ряда наборов дискретных направлений;
Разработанные алгоритмы реализованы в виде программного кода на языке FORTRAN.
Исследование характеристик поля излучения (радиационного потока, плотности лучистой энергии, дивергенции вектора плотности потока радиационной энергии) неоднородных плазменных объемов в зависимости от пространственных и угловых сеток, а так же от числа спектральных групп в оптическом диапазоне.
Тестирование разработанных методов на примерах сравнения результатов численных расчетов с имеющимися экспериментальными данными и с результатами, полученными другими методами решения уравнения переноса излучения.
Научная новизна.
Разработана методика применения метода дискретных ординат к расчету переноса теплового излучения потока вязкого, теплопроводящего газа через локализованную область плазмы в лазерном плазменном генераторе. Модель основана на уравнении переноса излучения в многогрупповом спектральном приближении. В качестве плазмообразующего газа исследован воздух и водород при атмосферном давлении. Проведен численный расчет радиационного нагрева внутренней поверхности лазерного плазменного генератора. Групповые и интегральные радиационные тепловые потоки на внутреннюю поверхность цилиндрического лазерного плазменного генератора были вычислены при помощи метода дискретных ординат на ортогональных структурированных сетках.
Произведено сравнение численных результатов метода дискретных ординат с методом дискретных направлений.
Разработан метод и представлены результаты численного моделирования радиационного нагрева задней поверхности космического аппарата MSRO (Mars Sample Return Orbiter) Европейского космического агентства. Для определения радиационных тепловых потоков разработан метод дискретных ординат на неструктурированных тетраэдальных сетках. Радиационная модель основана на уравнении переноса излучения в многогрупповом приближении.
Численный расчет выполнен для наиболее теплонапряженной точки предполагаемой траектории входа космического аппарата типа MSRO в атмосферу Марса содержащей 97% СОг и 3% N2 (массовые доли). Численные результаты получены для различных пространственных и угловых сеток. Произведено сопоставление результатов численных расчетов на структурированных и неструктурированных сетках. Предсказан уровень радиационных тепловых потоков к задней поверхности космического аппарата MSRO.
Для применения метода дискретных ординат для расчета поля излучения в объемах содержащих локализованные источники излучающего газа были вычислены веса SN угловых квадратур высоких порядков, вплоть до п=16. Практическая ценность работы.
Созданная расчетно-теоретическая модель, позволяет предсказывать характеристики радиационных полей в камерах лазерных плазменных генераторов, а таюке выполнять расчеты потерь радиационной энергии из воздушной и водородной лазерной плазмы энергетических устройств типа плазменных генераторов.
Вычисленные радиационные потоки на поверхность MSRO (Mars Sample Return Orbiter) входящего в атмосферу Марса, могут быть использованы для расчета тепловой защиты космического аппарата.
Разработанный расчетно-теоретический материал позволяет производить расчеты характеристик поля излучения в объемах произвольной формы, как для двумерных, так и для трехмерных геометрий. Созданный компьютерный код, основанный на методе дискретных ординат на треугольных (тетраэдальных) сетках, может быть использован для расчета поля излучения в радиационных газодинамических моделях в качестве дополнительного блока. Защищаемые положения.
Метод расчета переноса теплового излучения методом дискретных ординат на структурированных и неструктурированных сетках для двумерной, цилиндрической и трехмерной геометрий, особенностью которого является учет сильной неоднородности излучающего объема.
Квадратурные формулы высоких порядков точности для расчета теплового переноса излучения методом дискретных ординат в объемах с сильной локализацией плазмы.
Результаты расчета спектральных и интегральных характеристик поля излучения низкотемпературной плазмы в воздушном и водородном лазерном плазменном генераторе методом дискретных ординат.
Результаты расчета тепловых радиационных потоков к задней поверхности спускаемого космического аппарата при входе в атмосферу Марса.
Апробация работы.
Основные результаты диссертационной работы докладывались на следующих общероссийских и международных конференциях: XLIV научная конференция МФТИ, "Современные проблемы фундаментальных и прикладных наук", Москва, 23-30 ноября 2001
3 российская национальная конференция по теплообмену, Москва, 21-25 октября 2002 XLV научная конференция МФТИ, "Современные проблемы фундаментальных и прикладных наук", Москва, 29-30 ноября 2002
34th AIAA Plasmadynamics and Lasers Conference, Orlando, Florida, 23-26 June 2003
International Workshop on Radiation of High Temperature Gases in Atmospheric Entry, Lisbon, Portugal, 8-10 October 2003
IV международный симпозиум по радиационной плазмодинамике, Москва, 22-24 октября 2003 XL VI научная конференция МФТИ, "Современные проблемы фундаментальных и прикладных наук", Москва, 28-29 ноября 2003
40th AIAA Aerospace Sciences Meeting & Exhibit, Reno, NV, 8-11 January 2004 XLVII научная конференция МФТИ, "Современные проблемы фундаментальных и прикладных наук", Москва, 29 ноября 2004
38 AIAA Thermophysics Conference, Toronto, Ontario, Canada, 6-9 June 2005
Материалы, изложенные в диссертации, опубликованы в работах:
Филипский М.В., Суржиков СТ. Решение двумерной стационарной спектральной задачи переноса теплового излучения в неоднородной низкотемпературной плазме методом конечного объема // Труды XLIV научной конференции МФТИ, Изд-во МФТИ, 2001 г. С.47.
Филипский М.В. Решение двумерной спектральной задачи переноса теплового излучения в неоднородной низкотемпературной плазме методом конечного объема // Труды третьей Российской национальной конференции по теплообмену, Москва, Изд-во МЭИ, 2002 г. С.142-143.
Филипский М.В., Суржиков СТ. Метод конечного объема применительно к различным видам геометрии // Труды XLV научной конференции МФТИ, Изд-во МФТИ, 2002 г. С.12. Filipskiy M., Mokrov M., Surzhikov S., Capitelli M., Colonna G. Radiative Heating of Internal Surface of Hydrogen Laser Supported Plasma Generator II AIAA Paper 2003-4037, 34th Plasmadynamics and Lasers Conference, 23-26 June 2003, Orlando, Florida, P.I 1. Filipskiy M., Mokrov M., Surzhikov S., Capitelli M., Colonna G. Prediction of Radiative Heating of Internal Surfaces of Hydrogen and Air Laser Plasma Generators Intended for Aerospace Applications II 1 International Workshop on Radiation of High Temperature Gases in Atmospheric Entry; 8-10 October 2003, Lisbon, Portugal (ESA SP-533, December 2003), P.l 1-18.
Филипский M.B., Суржиков СТ. Решение двумерной осесимметричной задачи переноса теплового излучения в неоднородном объеме методом дискретных ординат // VI международный симпозиум по радиационной плазмодинамике, сборник научных трудов, М.: НИЦ "Инженер", 2003 г. С. 125-127.
Филипский М.В., Суржиков СТ. Применение метода конечного объема для решения двумерной цилиндрической задачи переноса теплового излучения в лазерно-плазменном ускорителе // VI международный симпозиум по радиационной плазмодинамике, сборник научных трудов, М.: НИЦ "Инженер", 2003 г. СЛ28-130.
Филипский М.В. Нахождение радиационных потоков на стенке воздушного плазменного генератора методом дискретных ординат // Труды XLVI научной конференции МФТИ, Изд-во МФТИ, 2003 г. С24. Filipskii М, Surzhikov S. Numerical Simulation of Radiation Heat Transfer in Plasma Generation II AIAA Paper 2004-0988, 40th Aerospace Sciences Meeting and Exhibit, 8-11 January 2004, Reno, NV, P. 12.
Филипский М.В. Модифицированный метод дискретных ординат применительно к криволинейной сложной геометрии // Труды XLVII научной конференции МФТИ, Изд-во МФТИ, 2004 г. С.31-32.
И.Филипский М.В. Расчет интегральных радиационных потоков на заднюю поверхность космического аппарата методом дискретных ординат // Труды XLVII научной конференции МФТИ, Изд-во МФТИ, 2004 г. С.33-34.
I2.Filipskiy М., Surzhikov S,, Discrete Ordinates Method for Prediction of Radiative Heating of Space Vehicales // AIAA Paper 2005-4948, 38th Thermophysics Conference, 6-9 June 2005, Toronto, Ontario, Canada, P.7. ІЗ.Филипский МБ., Суржиков СТ. Радиационный нагрев внутренней поверхности водородного и воздушного плазменного генератора // Вестник МГТУ им. Н.Э. Баумана. Сер. "Машиностроение". 2005.№2. С 3-19.
14. Филипский М.В., Суржиков СТ. Метод дискретных ординат для решения задач теплообмена излучением в сильно неоднородных средах // Препринт ИПМех РАН №781, 2005.
15.Филипский М.В., Суржиков СТ. Расчет радиационных потоков к поверхности космического аппарата с помощью метода дискретных ординат // Инженерно-физический журнал, том 79,2006.
Структура и объем работы.
Диссертация состоит из введения, девяти глав, заключения и списка литературы. Общий объем диссертации - 140 страниц, включая 95 рисунков и 4 таблицы. Список литературы содержит 57 наименований. Основное содержание работы.
Во введении обосновывается актуальность выбранной темы. Формулируются основные цели и содержание поставленных задач. Отмечена научная новизна и практическая значимость полученных результатов. Кратко сформулирована структура и содержание диссертационной работы.
Первая глава содержит литературный обзор по методу дискретных ординат. Приведена история развития метода и его вариаций: метод характеристик, метод конечного объема, Отмечено, что наиболее часто практический интерес представляет расчет переноса теплового излучения плазмы занимающую локализованную в объеме область, в то время как практически все расчетные методы теории переноса теплового излучения ориентированы на более или менее равномерное распределение источников излучения в исследуемой области.
Вторая глава посвящена применению метода дискретных ординат для решения уравнения переноса на ортогональных сетках.
В разделе 2.1 дана формулировка уравнения переноса теплового излучения в двух- и трехмерной постановках применительно к прямоугольной декартовой, цилиндрической и сферической систем координат.
В разделе 2.2 записаны конечно-разностные расчетные формулы интегрирования уравнения переноса теплового излучения в объемах разной геометрии с использованием ортогональных структурированных сеток. А так же приведены алгоритмы отыскания решения с использованием конечно-разностных уравнений. Особое внимание уделяется решению задач переноса в осесимметричной постановке.
Третья глава посвящена выбору значений весовых функций в уравнениях связи средних интенсивностей на граях со средней интенсивностью в центре ячейки для обеспечения положительности конечно-разностной схемы.
Метод дискретных ординат на ортогональных сетках, основные расчетные алгоритмы
В [13] применен SN- метод дискретных ординат для излучающей, поглощающей и анизотропно-рассеивающей среды заключенной в цилиндрическую полость с диффузно отражающими стенками. Функция рассеяния разлагалась в ряд по полиномам Лежандра, все расчеты производились с использованием S14 угловой квадратуры. Применение квадратур такого высокого порядка точности оправдано при расчете сред с высокой степенью анизотропии рассеяния. В работе исследуется влияние функции рассеяния на поток излучения падающего на стенку ограничивающей объем области.
В работе [14] использована стандартная формулировка SN- метода дискретных ординат для описания переноса излучения плазмы дугового разряда в аргоне при атмосферном давлении. Распределение поля температур в дуге считается известным и взято из ранее опубликованных работ. Спектральный коэффициент поглощения аргона рассчитывался по специализированному коду NASA. Для определения температурной зависимости спектрального коэффициента поглощения в условиях резкой температурной использовалась экспоненциальная интерполяция. Расчет производился в многогрупповом приближении. В работе исследовалась зависимость характеристик поля излучения плазмы от пространственных и угловых сеток, от температуры и количества спектральных групп.
В работе [15] произведено сравнение радиационных потоков в цилиндрических топках для различных угловых аппроксимаций низких порядков с точным решением. Сравнение показало, что полностью симметричное угловое распределение дает более точное решение, чем квадратура, основанная на кусочно-постоянном угловом распределении. Также произведено исследование чувствительности решения к изменению коэффициента поглощения.
В основе работы [16] лежит вывод уравнения переноса излучения в формулировке метода дискретных ординат применительно к криволинейной ортогональной системе координат. Решена задача сложного теплообмена и найдена температура в среде заключенной между двумя бесконечными эллиптическими цилиндрами. В [17] применен аналогичный подход для определения температурного поля в изогнутых симметричных пластинах. Результаты получены для эллиптических, гиперболических и параболических профилей.
Как уже отмечалось, успешным развитием метода дискретных ординат оказался метод характеристик. В методе характеристик численное интегрирование производится вдоль направления полета фотонов. В простейшем случае плоского слоя интегрирование вдоль характеристики проводится с использованием формального решения уравнения переноса. В прямоугольной декартовой системе координат характеристики являются прямыми линиями. Сложнее использовать метод характеристик в криволинейных системах координат. Владимировым B.C. был предложен алгоритм, в котором уравнение переноса в одномерной сферической геометрии интегрируется вдоль направления характеристики. Причем все характеристики проводятся параллельно некоторому направлению. Число угловых точек на каждом радиусе определяется числом пересечения окружности этого радиуса с характеристиками. Решение, полученное таким способом, сходится к точному решению уравнения переноса при уменьшении шага расчетной сетки по радиусу. Каждая характеристика может рассчитываться независимо от других, что делает данный алгоритм удобным для расчета на многопроцессорных ЭВМ. Однако наличие жесткой связи между пространственной и угловой сеткой делает алгоритм Владимирова B.C. не экономичным. Часто для численного решения задач переноса излучения требуется значительно менее подробная сетка по углу, чем по пространственной переменной. Геометрия сетки такова, что на радиусах, удаленных от центра шара, располагается значительное число угловых точек сетки, а вблизи центра их число мало.
В многомерных задачах узлы пространственной сетки могут не пересекаться с характеристиками. В методе характеристик с интерполяцией для нахождения интенсивности в узле ячейки уравнение переноса интегрируют вдоль характеристики от грани ячейки до узла, в котором ищется интенсивность. Интенсивность на грани ячейки отыскивается интерполяцией через известные значения интенсивностей в узлах ячейки принадлежащих грани. В работе [18] предложен один из методов интерполяции, позволяющий упростить логическую схему счета. Для криволинейных ячеек производится спрямление границ. Приведены расчетные формулы для разных типов геометрии. В сравнении со стандартным SN- методом метод характеристик с интерполяцией потребляет времени на 10-20% больше, зато имеет запас точности на более грубых сетках и не дает осциллирующего решения, чего нельзя сказать об исходном SN- методе.
Для двумерных областей в работе [19] применен метод характеристик с интерполяцией. Получены соотношения в предположении, что функция источника линейно изменяется вдоль характеристики. В работе применена линейная и кубическая интерполяция. Для параллелепипеда произведенное сравнение со стандартным методом дискретных ординат и с методом конечного объема показало более высокую точность метода характеристик. Причем при использовании кубической интерполяции ошибка расчета ниже, чем у линейной. В этой же работе выполнены расчеты для сред изотропно-рассеивающих и без рассеяния для разных оптических длин. Данный метод позволяет вести расчет на криволинейных сетках. В работе [20] использовался SN- метод для исследования переноса излучения в нерассеивающей среде (пары воды) в одномерной геометрии, ограниченной абсолютно черными стенками. Одномерная формулировка задачи позволила использовать метод характеристик для нахождения интенсивностей на гранях ячеек. Вычислялся полный тепловой поток на стенку при разных профилях температуры внутри области. Представлены результаты расчетов и произведено сравнение с результатами полученными методом потоков. Работа [21] является логическим продолжением работы [20], выполнена теми же авторами. Разработанный подход распространен на излучающую и поглощающую среду ограниченную диффузно-отражающими стенками. В указанной работе установлено, что для используемых параметров среды (пары воды) для достижения высокой точности решения хватает двукратного учета отражения на стенке. В случае оптически тонких сред необходимо учитывать большее число отражений от стенок.
Интегрирование уравнения переноса излучения в формулировке метода дискретных ординат позволяет связать среднюю интенсивность по объему ячейки со средними интенсивностями на гранях. В SN- методе указанные уравнения связи выбираются из условия непрерывности решения в ячейки. В работе [22] применялся метод характеристик для нахождения неизвестных интенсивностей на гранях через известные. В расчетной двумерной области вводится произвольная треугольная сетка. Внутри каждой ячейки характеристики среды считаются постоянными. В работе рассматривается задача переноса излучения с теплопроводностью, то есть задача сложного теплообмена для двумерных областей. Рассеивающая среда ограничена равномерно нагретыми диффузно отражающими стенками. Приведены радиационные потоки на поверхностях и температурные профили внутри параллелепипедов, круглых и эллиптических бесконечных цилиндров.
Цилиндрическая геометрия
В работе [20] использовался SN- метод для исследования переноса излучения в нерассеивающей среде (пары воды) в одномерной геометрии, ограниченной абсолютно черными стенками. Одномерная формулировка задачи позволила использовать метод характеристик для нахождения интенсивностей на гранях ячеек. Вычислялся полный тепловой поток на стенку при разных профилях температуры внутри области. Представлены результаты расчетов и произведено сравнение с результатами полученными методом потоков. Работа [21] является логическим продолжением работы [20], выполнена теми же авторами. Разработанный подход распространен на излучающую и поглощающую среду ограниченную диффузно-отражающими стенками. В указанной работе установлено, что для используемых параметров среды (пары воды) для достижения высокой точности решения хватает двукратного учета отражения на стенке. В случае оптически тонких сред необходимо учитывать большее число отражений от стенок.
Интегрирование уравнения переноса излучения в формулировке метода дискретных ординат позволяет связать среднюю интенсивность по объему ячейки со средними интенсивностями на гранях. В SN- методе указанные уравнения связи выбираются из условия непрерывности решения в ячейки. В работе [22] применялся метод характеристик для нахождения неизвестных интенсивностей на гранях через известные. В расчетной двумерной области вводится произвольная треугольная сетка. Внутри каждой ячейки характеристики среды считаются постоянными. В работе рассматривается задача переноса излучения с теплопроводностью, то есть задача сложного теплообмена для двумерных областей. Рассеивающая среда ограничена равномерно нагретыми диффузно отражающими стенками. Приведены радиационные потоки на поверхностях и температурные профили внутри параллелепипедов, круглых и эллиптических бесконечных цилиндров. Данный метод не имеет ограничений на форму области, так как все расчетные формулы записаны в декартовой системе координат. Авторы [22] утверждают, что метод не дает осцилляции решений.
В последующей работе тех же авторов [23] применен метод [7], позволяющий заметно повысить точность расчета сред с рассеянием, ограниченных локально нагретыми стенками. Если в работе [7] применялась двумерная ортогональная сетка, то в работе [23] уменьшение лучевого эффекта стало возможным в двумерных областях произвольной формы. Также в работе [24] приведены соотношения для трехмерной геометрии, используя подход развитый ранее в [22]. При численном сравнении с расчетными данными, даваемыми зональным методом для параллелепипеда содержащего нерассеивающую среду, метод с использованием характеристик дает более точное решение, чем стандартный метод дискретных ординат [9]. В работе приведены численные результаты его применения для расчета радиационных полей сред ограниченных стенками произвольной геометрии.
В работе [25] предложена схема численного решения задач переноса излучения в произвольных осесимметричных областях, которые могут содержать внутренние границы, разделяющие среды с различными показателями преломления, а также зеркально отражающие границы. Область разбивается на трехмерные пространственные ячейки, и уравнение переноса излучения интегрируется вдоль характеристик от одной грани ячейки до другой. Решение ищется во всей трехмерной области для лучей параллельных плоскости, которая содержит ось симметрии исследуемого объема и одну из главных координатных осей. Приводятся результаты для нерассеивающих сред, однако, судя по всему, предлагаемая схема достаточно легко обобщается на случай изотропного и анизотропного рассеяния.
Низкая точность стандартного SN- метода дискретных ординат для расчета поля излучения в областях, содержащих вакуумные полости, выдвигает необходимость разработки специальных методов нахождения решения. Тот факт, что в вакууме уравнение переноса излучения имеет аналитическое решение, позволил разработать специализированный VAC алгоритм [26]. В вакууме излучение распространяется без поглощения вдоль характеристик. Алгоритм является локальным, то есть при расчете каждой расчетной ячейки производится геометрический анализ только соседних ячеек. Благодаря своей локальности алгоритм VAC может рассчитывать перенос излучения без поглощения в объемах произвольной формы и легко включается в цикл расчета по разностной сетке. Таким образом, область с пустотами рассчитывается следующим образом: плотная среда рассчитывается каким-либо вариантом метода дискретных ординат, а вакуумные полости VAC алгоритмом.
В работах [27]-[29] приведен метод выбора весовых коэффициентов с помощью метода характеристик для SN- метода дискретных ординат. В [27], [29] исследованы балансные уравнения метода дискретных ординат на предмет возникновения нефизических осцилляции интенсивности. С помощью метода характеристик находится решение в ячейке [28]. Это решение используется для нахождения весов в балансовом уравнении переноса излучения метода дискретных ординат. В [28] приведены характеристические поверхности для расчетной ячейки в цилиндрической системе координат.
Результаты тестирования разработанных программных кодов
Точное решение уравнения (2.1.4) применительно к широкому кругу задач радиационного теплопереноса является трудоемким, а порой даже невозможным. Поэтому для нахождения решения уравнения переноса используют различные приближения метода. К настоящему времени разработана серия таких методов [49], но здесь будет рассмотрен только метод решения уравнения переноса в формулировке дискретных ординат на ортогональных сетках. Идея этого метода заключается в выборе конечного числа дискретных направлений и установлении связи между интенсивностями излучения, рассчитанными в указанных направлениях. Все интегралы по угловым переменным разлагаются в квадратурные формулы по конечному числу выбранных направлений. Таким образом, интегро-дифференциальная краевая задача (2.1.4) преобразуется в систему уравнений в частных производных, зависящих только от пространственных переменных.
Разобьем область угловых направлений Q на конечное число единичных векторов Qm (где т- номер дискретного направления, общим числом М) и каждому направлению присвоим некоторое весовое значение сот. Для лучшего понимания, представим точки т на поверхности сферы единичного радиуса, каждой точке соответствует некоторая площадь поверхности сферы. Сумма всех весов равняется площади поверхности сферы, то есть о т =Аж. Направляющие векторы имеют координаты Qm(/;m,//„,„,), где Т]т, ря ,- направляющие косинусы вектора Qm. Способу выбора угловых направлений (используется также термин угловых квадратур) посвящен отдельный раздел данной работы (см. шестую главу). Тогда уравнение (2.1.4) можно записать в виде системы уравнений переноса излучения для каждого направления С1т температура границы области, Я - нормаль к границе в точке г при г є Г. Интеграл рассеяния в правой части уравнения переноса заменяется квадратурной суммой. Индикатриса рассеяния (в случае анизотропии среды) разлагается по полиномам Лежандра [9]. В дальнейшем будем рассматривать (2.2.1) применительно к изотропно-рассеивающей среде (Ф = 1). Введем пространственную сетку. Для этого всю расчетную область разобьем на прямоугольные ячейки рис.2.2.1. Центр ячейки имеет координаты г іх,, у, zk ), где индекс / є (і, IX} означает разбиение по оси х, по оси у, кє(і,К2)- по оси z. Координаты грани ячейки обозначаются смещением соответствующего индекса на половину. Таким образом, объем ячейки с координатами r\xny,,zk} равняется Индекс р будет в дальнейшем означать, что величина относится к центру ячейки. Будем считать коэффициенты поглощения и рассеяния постоянными внутри каждой ячейки и равными кр и ар соответственно. Трехмерная геометрия. В работах [4], [8], [9] описан метод дискретных ординат применительно к трехмерной декартовой системе координат. Проинтегрируем (2.2.1) для случая изотропно-рассеивающей среды по объему элементарной ячейки VIJk (Vp) с координатами (x„yJtzk}
Решение (2.2.7) для каждого направления Qffl отыскивается посредством итераций, поскольку граничные условия и аппроксимация интеграла рассеяния зависят от неизвестных интенсивностеи. Для примера покажем алгоритм вычисления на первом шаге итераций в одной спектральной группе. Для того, чтобы начать вычисления полагают, что стенки расчетной области являются абсолютно черными, а коэффициент рассеяния равен нулю. В направлении Qm с положительными направляющими косинусами вычисления начинаются с ячейки с координатами (i,j,k) = (1,1,1). Именно для этого элементарного объема все три интенсивности на опорных гранях ячейки известны из граничных условий. С помощью (2.2.7) вычисляется интенсивность в центре ячейке. Используя (2.2.5) вычисляются интенсивности на конечных гранях, эти конечные грани являются опорными для граничащих ячеек. Таким образом, для соседних граней становится известными все три величины необходимые для решения (2.2.7). Процесс продолжается, пока интенсивности в каждой ячейке во всей расчетной области ни будут найдены. Для другого дискретного направления выбирается ячейка, в которой все три интенсивности на опорных гранях известны из граничных условий. Первая итерация заканчивается, когда все направления будут просканированы. На втором итерационном шаге расчетный алгоритм аналогичен приведенному выше, отличаясь только тем, что в граничных условиях (2.2.4), в слагаемом дающим вклад от отражения, и в уравнении (2.2.7), в аппроксимации интеграла рассеяния, неизвестные интенсивности берутся с предыдущей итерации. Итерации прекращаются при достижении заданной точности решения по любому из наперед заданных критериев.
Результаты численного моделирования водородного ЛПГ
После описания алгоритма решения уравнения переноса излучения методом дискретных ординат, осталось определить используемые угловые ординаты и соответствующие им веса. Методы расчета угловых ординат обсуждаются в работах [41]-[47], где приведены квадратурные разложения. Напомним, что в уравнение переноса излучения входит интеграл рассеяния по телесному углу 4я, в граничные условия входит интеграл диффузного отражения стенки по половине телесного угла 2тг. Некоторые характеристики поля излучения находятся также путем интегрирования по сфере (например, поток лучистой энергии в направлении Clr, плотность лучистой энергии). Поэтому, естественным является использование единого набора направлений для аппроксимации интегралов следующих типов
То есть задача сводится к нахождению точек на сфере и присвоению им соответствующих весов, чтобы как можно точнее брались интегралы (6.1). Согласно методу дискретных ординат [45] разложим интенсивность в ряд по ортогональным сферическим функциям где JV- порядок аппроксимации интенсивности сферическими функциями; є ошибка аппроксимации. Подставим (6.2) в (6.1), тогда интегрирование интенсивности сведется к интегрированию сферических функций. Для того чтобы уменьшить погрешность расчета интегралов (6.1) потребуем точного интегрирования сферических функций до порядка аппроксимации JV включительно квадратурными рядами порядка М. То есть, для всех n N, 1 п потребуем точного выполнения равенств 2. Направления Qm должны не зависеть от выбора пространственного базиса, т.е. при любом вращении на 90 вокруг координатных осей квадратурные формулы не должны изменяться. 3. Все веса должны быть положительными. Обычно на сфере в первом октанте выбирают направление Qm(a (3 у) и при последовательном вращении получают точки в первом октанте (а Р у), (а у р), (р а у), (р у а), (у а р), (у р а). Веса этих направлений должны быть одинаковыми. Существует три способа задания симметричных направлений. Когда все три координаты различны, возможны шесть направлений с идентичными весами. Для случая, когда две координаты вектора совпадают, получается три направления с одинаковыми весами. И когда вектор имеет все три одинаковые координаты, возможен только один случай. В остальных октантах получают направления, используя первый октант, путем вращения сферы на 90 вокруг главных осей. Например, в работе [46] рассмотрена так называемая TN квадратура (рис.6.1). Она получается путем разбиения плоскости г] + ц + % = \ в первом октанте на одинаковые треугольники. В записи TN, N означает на сколько частей разбивается сторона главного треугольника. Центральные точки треугольников соответствуют дискретным направлениям. Проекции данных треугольников на сферу единичного радиуса образует TN квадратуру. Веса такой квадратуры находят, как площади проекций треугольников на сфере. В работе [4] предложен выбор дискретных направлений, которые лежат на широтах сферы. Этот способ выбора направлений называется SN квадратурой, где N означает количество проекций на ось. На рис.6.2 представлена S8 квадратура. В одном октанте всего четыре проекции направлений Пт на ось. В силу симметрии относительно вращения, всего восемь проекций. Координаты на поверхности сферы даны формулой [4] От выбора первого направляющего косинуса / зависит распределение направлений по сфере. При малых /if направления расположены ближе к осям. При больших значениях $ — , направляющие косинусы группируются блшке к направлению Для нахождения весов квадратуры используют систему уравнений (6.5Ь), В силу симметрии (6.5а) будет выполняться автоматически. Для записи системы уравнений (6.5Ь), выберем половину сферы 0. Из симметрии относительно вращения вокруг оси следует, что интегралы с нечетными степенями /,, /2при направляющих косинусах rj и ц будут равняться нулю. Значения интегралов полученных перестановкой степеней имеют одни и те же значения. Рассмотрим пример расчета весов. Так как мы используем симметричные направления относительно вращений вокруг координатных осей, интегралы можно брать в первом октанте. В S8 квадратуре в первом октанте всего три вектора направления: П,(а,,Д,/,), Q.2(a2,a2,y2), Q3(a3,a3,a3). Остальные направления получаются после вращения.
