Линейное уравнение Больцмана: приближение, методы численного решения прямых задач и задач оптимизации, обобщение тема диссертации и автореферата по ВАК РФ 05.13.18, доктор наук Руколайне Сергей Анатольевич

• Руколайне Сергей Анатольевич
• доктор наук
• 2019

Оглавление диссертации доктор наук Руколайне Сергей Анатольевич

1.3.2 Оптимизация формы области

1.4 Неклассический перенос

2 приближение к линейному уравнению Больцмана

2.1 Модели диффузии

2.1.1 Уравнение диффузии

2.1.2 Телеграфное уравнение

2.1.3 Уравнение типа Джеффриса

2.2 Приближения к линейному уравнению Больцмана в рамках метода сферических функций

2.2.1 Метод сферических функций

2.2.2 Диффузионное приближение

2.2.3 Рм приближения

2.2.4 Р1/3 приближение

2.2.5 приближения

2.2.6 приближение

2.3 Задача Коши для уравнения в трехмерном пространстве

2.3.1 Локальный источник малой продолжительности

2.3.2 Локальное начальное распределение, источник отсутствует

2.4 Резюме

2.Л Приложение

2.Л.1 Вывод решения задачи Коши (2.34), (2.37)

3 Методы численного решения задач переноса излучения в областях с

зеркальными границами

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

3.1.1 Краткий исторический экскурс

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

3.1.3 Кусочно-квазилинейная интерполяция, основанная на триангуляции 1-го типа

3.1.4 Кусочно-квазилинейная интерполяция, основанная на триангуляции 2-го типа

3.1.5 Дискретизация граничного условия

3.2 Уравнение переноса излучения и граничные условия в случае осевой симметрии

3.2.1 Уравнение переноса

3.2.2 Другой способ вывода уравнения переноса

3.2.3 Условие на непрозрачной диффузной границе

3.2.4 Условия на прозрачной диффузной границе раздела сред

3.2.5 Условие на непрозрачной зеркальной границе

3.2.6 Условия на прозрачной зеркальной границе раздела сред

3.3 Упрощенная численная схема решения осесимметричных задач переноса излучения

3.3.1 Разбиение области

3.3.2 Дискретизация уравнения переноса

3.3.3 Дискретизация граничных условий

3.3.4 Алгоритм численного решения

3.3.5 Тестовые задачи

3.4 Численная схема решения осесимметричных задач переноса излучения

3.4.1 Разбиение области

3.4.2 Дискретизация уравнения переноса и граничных условий

3.4.3 Алгоритм численного решения

3.4.4 Тестовые задачи

3.4.5 Модификация численной схемы

3.5 Применение при моделировании роста полупрозрачных кристаллов

3.Л Приложения

3.Л.1 Вычисление интегралов по сферическим треугольникам, образованным дугами большого круга

3.Л.2 Интегрирование кусочно-квазилинейных функций первого типа по

сферическим треугольникам

3.Л.3 Построение квазилинейных функций второго типа на сферических

треугольниках

4 Оптимизация граничных значений и формы области в задачах переноса излучения

4.1 Регуляризация обратных задач оптимизации граничных значений

4.1.1 Постановка обратной задачи оптимизации граничных значений

4.1.2 Используемые методы регуляризации

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

4.1.4 Градиент функционала

4.1.5 Ограничения на искомое решение

4.1.6 Численная реализация

4.1.7 Модельная задача

4.1.8 Тестовые задачи

4.1.9 Регулярное решение осесимметричных обратных задач оптимизации граничных значений

4.1.10 Резюме

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

4.2.1 Постановка обратной задачи геометрической оптимизации

4.2.2 Градиент целевого функционала

4.2.3 Вывод интегрального тождества, определяющего субстанциальную производную

4.2.4 Пример

4.2.5 Вывод интегрального тождества, определяющего вариационную производную

4.2.6 Вывод сопряженной задачи и вычисление градиента целевого функционала

4.2.7 Общая схема процедуры вычисления градиента целевого функционала

4.2.8 Градиент целевого функционала в случае, когда оптимизируемая поверхность задана конечным числом параметров

4.2.9 «Двухмерные» области, в которых оптимизируемая граничная поверхность представляет из себя многогранник

4.2.10 Модельная задача

4.2.11 Общая схема процедуры вычисления градиента целевого функционала (продолжение)

4.2.12 Тестовая задача с диффузными границами

4.2.13 Тестовые задачи с диффузно-зеркальными границами

4.2.14 Резюме

4.A Приложение

4.A.1 Некоторые формулы дифференциального исчисления на поверхностях

5 Обобщенное линейное уравнение Больцмана, описывающее неклассический перенос, и асимптотическое приближение к нему

5.1 Неклассический перенос частиц: модель

5.2 Обобщенное линейное уравнение Больцмана

5.3 Интегральные уравнения для плотностей

5.4 Асимптотическое решение при малой средней длине пробега

5.4.1 Предположения

5.4.2 Уравнение и представление решения

5.4.3 Внешнее решение

5.4.4 Внутреннее решение

5.4.5 Начальные условия

5.4.6 Резюме

5.5 Диффузионное приближение к односкоростному обобщенному линейному уравнению Больцмана

5.A Приложения

5.A.1 Асимптотическое разложение ядра ф

5.A.2 Разрешимость интегральных уравнений Фредгольма второго рода (5.27)

Заключение

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

Введение диссертации (часть автореферата) на тему «Линейное уравнение Больцмана: приближение, методы численного решения прямых задач и задач оптимизации, обобщение»

Введение

Актуальность. Линейное кинетическое уравнение Больцмана, называемое также нестационарным уравнением переноса излучения или уравнением переноса нейтронов, описывает, например, перенос теплового излучения (радиационный теплоперенос — перенос тепловой энергии электромагнитными волнами) [32, 37, 50, 143, 196, 211] и перенос нейтронов [16, 18, 29, 51, 101, 177]. Кроме того, линейное уравнение Больцмана используется для построения моделей переноса или передвижения биологических объектов или служит основой для построения таких моделей [61-63, 70, 71, 103, 136, 205-207, 277]. Линейное уравнение Больцмана — это интегро-дифференциальное уравнение с частными производными относительно плотности частиц в фазовом пространстве, которая в общем случае зависит от семи независимых переменных: время, три пространственные координаты, две угловые координаты (направление движения) и энергия (или скорость) частиц.

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

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

говоря, не распространяются на области более сложной формы.

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

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

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

Научная новизна:

• приближение к линейному уравнению Больцмана предложено в качестве уточнения диффузионного приближения.

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

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

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

• Предложено обобщение линейного уравнения Больцмана на случай произвольного распределения длины свободного пробега частиц.

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

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

Положения, выносимые на защиту:

1. Dl приближение предложено в качестве уточнения диффузионного приближения к линейному уравнению Больцмана.

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

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

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

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

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

Апробация результатов. Основные результаты диссертационной работы докладывались на Девятой Российской национальной конференции по росту кристаллов (Москва, 2000), The Third International Symposium on Radiative Transfer (Antalya, Turkey, 2001), The Fourth International Conference on Single Crystal Growth and Heat & Mass Transfer (Обнинск, Россия, 2001), на Десятой Российской национальной конференции по росту кристаллов (Москва, 2002), на Третьей Российской национальной конференции по теплообмену (Москва, 2002), Eurotherm Seminar «Computational Thermal Radiation in Participating Media» (Mons, Belgium, 2003), The Fourth International Symposium on Radiative Transfer (Istanbul, Turkey, 2004), Eurotherm Seminar «Computational Thermal Radiation in Participating Media II» (Poitiers, France, 2006), на Четвертой Российской нацио-

нальной конференции по теплообмену (Москва, 2006), на Пятой международной конференции по обратным задачам (Казань-Москва, 2007), The Fifth International Symposium on Radiative Transfer (Bodrum, Turkey, 2007), на Международной конференции «Days on Diffraction 2013» (С.-Петербург, 2013), на Международном симпозиуме «Systems Biology and Bioinformatics» (С.-Петербург, 2016), на семинарах в Физико-техническом институте им. А. Ф. Иоффе и Санкт-Петербургском политехническом университете.

Публикации. По теме диссертации опубликовано 27 статей, из них 16 статей в рецензируемых научных журналах, индексируемых Scopus и/или Web of Science и/или входящих в перечень ВАК [10, 27, 34, 35, 76, 107, 226, 228, 230, 233, 235, 239, 241, 243, 244, 281], 10 статей в трудах научных конференций и научных сборниках [33, 73, 227, 229, 231, 232, 234, 236-238] и 1 препринт [242].

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

Благодарности

Выражаю глубокую благодарность В. С. Юфереву, привлекшему меня к решению

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

Выражаю искреннюю признательность О. Н. Буденковой, М. Г. Васильеву, В. М. Мамедову и О. И. Чистяковой, без сотрудничества с которыми настоящая диссертация была бы невозможна.

Глава 1

Линейное уравнение Больцмана и связанные с ним проблемы

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

1.1 Приближения к линейному уравнению Больцмана

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

В сильнорассеивающих и слабопоглощающих средах, то есть средах с малым числом Кнудсена (отношением длины свободного пробега к характерному пространственному масштабу) и доминирующим рассеянием, широко используется диффузионное приближение к линейному уравнению Больцмана [101, 167, 196, 211]. Однако оно слиш-

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

Для построения приближений к линейному уравнению Больцмана широко используется метод сферических функций [121]. Решение уравнения Больцмана (плотность частиц в фазовом пространстве) представляется в виде ряда по сферическим функциям, в результате получается бесконечная система уравнений в частных производных относительно коэффициентов разложения (моментов): система моментных уравнений. Для того, чтобы приблизить решение уравнения Больцмана своими моментами до какого-то порядка, требуется так называемая процедура замыкания, поскольку каждое уравнение из бесконечной системы связывает моменты разных порядков. Существуют различные процедуры замыкания, исходящие из различных принципов, см. [121], которые приводят к различным приближениям к уравнению Больцмана. В простейшем и наиболее используемом случае бесконечная система просто «обрезается», то есть учитываются уравнения до заданного порядка, включительно, а моменты более высоких порядков, входящие в эти уравнения отбрасываются. В результате получаются Рм приближения. Наиболее используемое из Рм приближений — Р1 приближение.

Р1 приближение [101, 189, 211] «исправляет» диффузионное: оно описывается телеграфным уравнением (волновым уравнением с затуханием), которое относится к гиперболическому типу. Поэтому скорость распространения частиц в Р1 приближении конечна. Заметим, что телеграфное уравнение было предложено в качестве замены не только уравнению диффузии, но и уравнению теплопроводности [30, 150-152, 274]. Однако принципиальный недостаток Р1 приближения заключается в том, что оно дает неверную скорость распространения частиц (в трехмерном пространстве), а именно ь/л/З, где V — истинная скорость частиц. Простая модификация Р1 приближения, называемая Р1/3 приближением, была предложена в статье [204]. Р1/3 приближение также описывается телеграфным уравнением, но с одним отличающимся коэффициентом. Р1/3 приближение не только правильно описывает скорость распространения частиц, но и во многих случаях лучше, чем Р1 приближение [196, 198, 204], см. также [203]. Однако как Р1/3, так и Р1 приближение имеют один общий принципиальный недостаток: для телеграфного уравнения не выполняется принцип максимума, то есть телеграфное уравнение не сохраняет неотрицательность решений, и поэтому решения задач для него могут быть нефизичны [162, 212]. То же самое имеет место и для асимптотических Р1 приближений, описываемых телеграфным уравнением [133, 134].

Рм приближения при N > 1 имеют тот же принципиальный недостаток, что и Р1 приближение: они не сохраняют неотрицательность решений (имеют вид распро-

страняющихся волн), см. [121]. Поэтому PN приближения не могут служить заменой диффузионному приближению.

Существуют также упрощенные PN (simplified PN, SPN) приближения [67, 122124, 188, 196], но и они имеют свои ограничения.

Недавно были предложены DN, N ^ 1, приближения к односкоростному уравнению Больцмана [248], причем диффузионное приближение можно рассматривать как D0 приближение. Заметим, что модели замыкания моментов уравнения Больцмана (в частности Dn модель) предназначены, вообще говоря, для того, чтобы аппроксимировать само уравнение Больцмана с меньшими вычислительными затратами, (а не получать приближения к нему для сильнорассеивающих сред), см., например, [121]. В статье [248] DN решения сравнивались численно с PN решениями в двух тестовых задачах: в одномерном и двухмерном пространствах. Сравнение показало, что «DN модель дает значительно лучшие результаты при меньших N, чем PN модель» [248]. В статье [121] различные аппроксимации к уравнению Больцмана, в частности PN и DN приближения, сравнивались в двухмерной эталонной задаче (линейный источник излучения). Было обнаружено, что качественное поведение DN и PN_1 решений похоже, но Dn решения ослабляют в некоторой степени колебательный характер Р^-1 решений. Таким образом, двойником Р1 приближения является D2 приближение, а двойник у D1 приближения отсутствует, поскольку не существует Р0 приближения.

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

1.2 Методы решения задач переноса излучения 1.2.1 Уравнение переноса излучения

Нестационарное уравнение переноса излучения это линейное уравнение Больцмана, описывающее перенос теплового излучения (перенос тепловой энергии электромагнитными волнами) [32, 37, 50, 143, 196, 269] (см. также [193, 194]). Оно имеет вид

1 д1

+ П •УI, + /Зи 1и = ^ 5 4 + к„ 1ъ>и, (1.1)

С ОЪ '

излучения, г — пространственные переменные,

О = (sin в cos sin в sin cos в)

— направление излучения, t — временная переменная, v — частота, играющая роль энергии (и = с/Х, X — длина волны, энергия кванта электромагнитной волны (фотона) равна hu, где h — постоянная Планка), ки — коэффициент поглощения, asu — коэффициент рассеяния, = kv + as v — коэффициент ослабления,

5Iv =SIv(г, П) = ± Í Фу(ü, ü')Iu(r, ü')dü'

J §2

— интеграл рассеяния (d.O = sin 0 &в d^), Фv(Q, Q') — фазовая функция (или индикатриса) рассеяния, для которой выполняется нормировка

1

— Í Фу (П, ü') dQ = 1

J §2

(то есть

-1 Ф„(П, О'

— вероятность того, что излучение, имеющее частоту V, и распространяющееся в направлении О!, будет рассеяно внутри элементарного телесного угла в направлении

П),

I =1 т-- 2Нр3

1Ъ, V — 1Ъ, V ) —

С2

exp ( S - 1

— монохроматическая интенсивность излучения черного тела (функция Планка) [32, 196], Т — Т(г) — температура среды, Н — постоянная Планка, кв — постоянная Больц-мана.

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

П •У/, + & 1и — а8 ,, 51и + ки 1ъ ,,, (1.2)

где 1и — 1и(г, О) — стационарная интенсивность излучения.

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

В уравнение переноса (1.2) (так же, как и в нестационарное уравнение (1.1)) ча-

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

П •У! + ¡31 = а^ I + к1ъ, (1.3)

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

1.2.2 Методы численного решения задач переноса излучения

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

К стохастическим методам относится метод Монте Карло [140, 141, 143, 196], который в применении к задачам переноса излучения основан на стохастическом моделировании траекторий фотонов. Если набрана достаточно большая статистика (смоделировано большое количество траекторий), то метод позволяет получать решения с высокой степенью статистической точности независимо от сложности задачи переноса излучения, что делает его универсальным средством для получения эталонных решений. Но стохастическая природа метода является и его недостатком, поскольку получение решений с высокой точностью требует больших вычислительных затрат. Это делает метод непригодным для решения задач, в которых перенос тепла излучением является лишь одним из видов тепло- и массопереноса.

К детерминистическим методам относятся метод сферических функций (Рм приближение) [101, 143, 192, 196, 211, 269], метод дискретных ординат приближение) [50, 88, 143, 196, 269] и метод конечных объемов [85, 86, 88, 143, 196, 216], зональный метод [143, 185, 196], метод дискретного переноса [89, 143, 183, 196, 252], и др.

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

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

Метод дискретных ординат (метод «многих потоков») был предложен Чандра-секаром [50] как обобщение так называемого метода «двух потоков» Шустера [249] и Шварцшильда [250] и сначала нашел применение в теории переноса нейтронов [80, 82, 168, 172, 174, 177], а затем при решении задач переноса теплового излучения [110— 112, 114, 145, 217, 218, 245-247, 262, 263].

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

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

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

Метод дискретного переноса комбинирует особенности метода дискретных ординат, зонального метода и метода Монте Карло.

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

Одним из наиболее удачных способов пространственной дискретизации в методе дискретных ординат является метод характеристик. Для решения задач переноса излучения метод характеристик был впервые предложен в статье [11]. Первоначально он развивался для решения задач переноса нейтронов [29, 168], где рассматривались задачи в декартовой геометрии [53, 69, 186] или в упрощенных сферической и цилиндрической геометриях [29]. Следует отметить, что задачи преноса нейтронов существенно отличаются от задач переноса теплового излучения, поскольку в них отсутствует отражение и преломление на границах. В применении к задачам переноса излучения подход, практически идентичный методу характеристик, был предложен в статьях [246, 247], в которых рассматривались задачи в двухмерной и трехмерной декартовых геометриях, соответственно. В декартовых координатах метод характеристик концептуально прост, поскольку дифференциальный оператор в уравнении переноса действует только на пространственные переменные (угловые переменные являются параметрами), и поэтому характеристики это прямые линии.

Список литературы диссертационного исследования доктор наук Руколайне Сергей Анатольевич, 2019 год

с источником

где

р(r, Pг(r)X[0,в](t),

ХадСО

-, 0 < в, 0, > .

Заметим, что начальный момент времени равен Ь = 0, источник равен нулю при Ь > в. Начальные условия для уравнения (2.48) нулевые, то есть

ди

Я=о = 0, Ж

0.

(2.49)

г=о

Замена Ь в задаче (2.48), (2.49) на Ь + 9 приведет к задаче (2.34), (2.37) Преобразование Фурье задачи (2.48), (2.49) дает

д2Ти

^ + [1+Т ^ +

дТи

+ [Р>! + ^2) \1\2 + 7]

Тр е[ (1 + ^^ХмЮ + тх'мМ

дТи

Ти\

г=о

0,

дЪ

0,

| е Е3, г> 0, (2.50)

(2.51)

г=о

где

Тр£ = Тр£ (|) = е

-в2|$|2/2

Преобразование Лапласа задачи (2.50), (2.51) дает

ТСи(£,з) = Тр£Сх[о,в] Т +,| 1 1Т(5 - Ау)(8 - Л2)

1 +8 + ^

5 — А, 5 — А,

_ ^р£СХ[о,в] [ Ау + 7 А2 + 7

Ат — А2

в — А2 в — А1

)

где

Си(■, в)

и(-,г) е-3* &

— преобразование Лапласа, и Ау 2 характеристические значения (2.39). Выполняя обратное преобразование Лапласа, получим

Т и(£,Ь)

Тр£

А- - А2

(А- + 7) е^

1 - е-л2б

- (А2 + 7) еЛх*

1 е-

Хув

г > в.

1

1

ос

0

Замена £ в правой части на £ + в дает формулу (2.38).

Глава 3

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

В разделе 3.1 построены квадратурные схемы, названные PQLA^ (Piecewise QuasiLinear Angular, кусочно-квазилинейные угловые), N = 1, 2,... , которые используют кусочно-квазилинейную интерполяцию на единичной сфере. Такие схемы дают квадратурные формулы на сфере для вычисления интеграла рассеяния и объемной плотности энергии излучения, и квадратурные формулы на полусферах для вычисления потока излучения на границах. Благодаря аналитическому представлению угловой зависимости, формулировка граничных условий в случае зеркального отражения и преломления на границе в рамках метода дискретных ординат становится простой. Кроме того, аналитическое представление угловой зависимости интенсивности излучения может быть использовано, если, например, требуется вычислять поток излучения в некотором пространственном угле.

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

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

В разделе 3.3 описана упрощенная численная схема решения задач переноса излу-

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

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

В разделе 3.5 описано применение разработанных численных схем при моделировании роста полупрозрачных кристаллов германата висмута из расплава низкоградиентным методом Чохральского, что позволило объяснить особенности роста этих кристаллов.

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

Этот раздел основан на результатах статьи [226].

3.1.1 Краткий исторический экскурс

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

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

На практике наиболее часто используются четыре типа квадратур: Ньютона-Котеса, Гаусса, Маркова и Чебышева. В квадратурах Ньютона-Котеса узловые точки заданы заранее, а квадратурные веса должны быть определены из условия, чтобы квадратура была точна для всех полиномов до некоторой степени. В квадратурах Гаусса узловые точки заранее не фиксированы, их координаты вместе с весами должны быть определены из такого же условия. В квадратурах Маркова заранее задана лишь часть узлов, а остальные узлы и все веса должны быть определены из условия. В квадратурах Чебышева все веса заданы и равны 1/Ы, где N — число узловых точек, а их координаты должны быть определены.

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

В случае метода дискретных ординат представляется естественным [80, 174] строить квадратурные схемы, инвариантные относительно группы симметрии октаэдра см. рис. 3.1. Эта группа оставляет инвариантным октаэдр, вписанный в единичную сферу, вершины которого лежат на координатных осях. Группа содержит (помимо тождественного преобразования) следующие преобразования: вращения на углы, кратные ^/2, вокруг каждой из координатных осей; вращения на углы, кратные 2^/3, вокруг осей, проходящих через центры граней октаэдра; вращения на углы, кратные п, вокруг осей, проходящих через середины ребер октаэдра; инверсию.

Рисунок 3.1. Октаэдр.

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

1. (±1,0, 0), (0, ±1, 0), (0,0, ±1) — вершины октаэдра;

2. (±1/л/3, ±1/л/3, ±1/л/3) — проекции центров граней октаэдра на сферу;

3. (±1/^2, ±1/^2, 0), (±1/^2, 0, ±1/^2), (0, ±1/^2, ±1/^2) — проекции середин ребер октаэдра на сферу;

4. (±Рк, ±Як, 0) (±Рк, 0, ±Як)> (0, ±Рк, ±Qk), где Рк + Як = 1, Рк = 0 Як = 0, Рк = Як, — центральные проекции этих точек лежат на ребрах октаэдра;

5. (± г^ ± г к, ± « к), (± гк, ± 5 к, ± гкX (± 5 к, ± гк, ± гкX где 2 гк + 5 к = 1 гк = 0, 5 к = 0

к = к , — центральные проекции этих точек лежат на высотах граней октаэдра;

6. (±а, ±@к, ±Ъ), (±ак, ±Ъ, ±Рк)> (±рк, ±ак, ±Ъ)> (±рк, ±1к, ±ак),

(± 1к,±а,±РкX (± 1к,±Рк,±а)> где ак + Рк + чк = 1, а = 0 Рк = 0 ъ = 0,

а к = Рк = 1к, — центральные проекции этих точек лежат на гранях октаэдра.

Узловые точки, принадлежащие разным классам, неэквивалентны, то есть не могут быть переведены друг в друга преобразованиями симметрии. Все узловые точки, принадлежащие каждому из 1-го, 2-го и 3-го класса, эквивалентны. Каждый из 4-го, 5-го и 6-го классов состоит из непересекающихся подмножеств эквивалентных узловых точек (например, узловые точки р1 и р2, р1 = р2, р\ + рк, = 1, образуют два неэквивалентных подмножества 4-го класса). Отметим, что множество узловых точек, образующих некоторую квадратуру, не обязано содержать точки из всех классов.

Одно из основных утверждений статьи [38] содержится в следующей теореме: Теорема. Пусть Ь — конечномерное пространство функций, инвариантное относительно преобразований группы С. Квадратура, инвариантная относительно группы С, точна для всех функций f Е Ь тогда и только тогда, когда она точна для всех функций f Е Ь, инвариантных относительно группы С.

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

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

Уп,т(#,<р) = РПН) (сое')ет = -п,...,п, п = 0,1, 2,...,

где Рп — присоединенные функции Лежандра, (',<р) — сферические координаты. Пусть Ь=п — пространство, образованное сферическими функциями Упт порядка п, т = - п, . . . , п, и

Ь<п = Ь=о и Ь=1 и... и Ь=п.

Размерность пространства Ь=п равна 2п +1, размерность пространства ¿<га равна (п + 1)2. Оба пространства, Ь=п и Ь<п, инвариантны относительно группы С\,2и. Сферические функции Упт, т = -п,...,п, являются линейными комбинациями полиномов П%ПЬУП*, а + Ь + с < п, где

üx = sin $ cos Qy = sin $ sin Qz = cos

поэтому пространство L<n может быть определено как пространство, образованное полиномами QaxQbyQcz, а + b + с < п. Таким образом число линейно независимых полиномов степени, не большей п, равно размерности пространства L<n, то есть (п + 1)2.

Согласно приведенной теореме квадратура на сфере, инвариантная относительно группы Gyjjj, точна для всех полиномов степени, не большей п, если она точна для всех таких полиномов, инвариантных относительно группы Gy2ii. Это утверждение значительно упрощает построение квадратур, поскольку число полиномов степени, не большей п, инвариантных относительно группы Gy 211, значительно меньше (п +1)2 [38]. Например, для п =11 число таких полиномов равно 7.

Используя эту теорему, В. И. Лебедев в статьях [23-25] построил квадратуры Маркова, инвариантные относительно группы Gyjjj и точные для полиномов степени, не большей п, для п = 9,11,..., 29. Он построил также квадратуры Чебышева для п = 11,15 [24].

Приведем в качестве примера квадратуры Маркова 11-го и 17-го порядков из [23] (здесь и ниже мы обозначаем их символами L11 и L17, соответственно), квадратурные веса, соответствующие узловым точкам п-го класса, обозначаются символами (равенство квадратурного веса нулю означает, что соответствующие узловые точки отсутствуют):

L11 (число узловых точек равно 50):

w(1) = 4/315 w 0.0127; w(2) = 27/1280 w 0.0211; w(3) = 64/2835 w 0.0226;

w(4) = 0; w55) = 114/725760 w 0.0202, r1 = 1/УТТ w 0.3015113; w(6) = 0.

L17 (число узловых точек равно 110):

w(1) w 0.0038282705; w(2) w 0.0097937375; w(3) = 0;

pi w 0.87815891, w{4) w 0.0096949964;

r1 w 0.18511564, w[5) w 0.0082117373;

г2 w 0.39568947, wf] w 0.0095954713;

r3 w 0.69042105, W¡5) w 0.0099428149; w(6) = 0.

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

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

Интересно сравнить квадратуры Лебедева с квадратурами, построенными Кохом и др. в статье [160], где был использован близкий подход: было предложено строить квадратуры Гаусса и Маркова, инвариантные относительно группы GyIH, точные для полиномов на главных полусферах (х > 0, х < 0, у > 0, у < 0, z > 0, z < 0). Были использованы узловые точки второго, пятого и шестого классов, то есть точки не лежащие на координатных плоскостях. Заметим, что практически все квадратуры, используемые при решении уравнения переноса излучения, используют узловые точки только из этих классов.

Таблица 3.1 сравнивает относительную точность квадратур Лебедева L11 и L17 с квадратурами DCT020 - 1246, DCT020 - 2468, DCT111 - 1246810 и DCT111 - 24681012 из статьи [160] (эти квадратуры обозначаются символами DCT020a, DCT0206, DCT111" и DCT1116, соответственно) при вычислении моментов f 0 ^а^у^С dQ на полусфере z > 0. Заметим, что все квадратуры, инвариантные относительно группы GyIH, точны для момента (0, 0,0) и моментов (0, 0, 2) на главных полусферах.

Таблица 3.1. Относительные ошибки при оценивании моментов fz>0 ^С на полусфере z > 0.

Квад- Число Относительные ошибки (в %) для моментов (а, Ь, с)

ратура узлов

(0, 0, 1) (0, 0, 3) (0, 2, 1) (0, 0, 4) (0, 2, 2) (0, 0, 5) (0, 4, 1) (0, 2, 3) (2, 2, 1)

DCT020a 48 0 0.28 -0.28 0 0 -0.28 -3.60 1.40 6.46

DCT0206 48 2.20 -0.22 4.60 0 0 0.05 7.04 -0.80 8.07

L11 50 -1.23 0 -2.48 0 0 0 -2.48 0 -7.50

DCT111a 80 0 0.02 -0.04 0 0 0 -1.44 -0.12 3.87

DCT1116 80 1.30 -0.08 2.64 0 0 0.02 3.92 -0.24 4.50

L17 110 -0.56 0 -1.12 0 0 0 -1.44 0.01 -2.49

Точные значения ж ж/2 ж/4 2ж/5 2ж/15 ж/3 ж/8 ж/12 ж/ 2

Подчеркнем еще раз, что квадратуры Ь11 и Ь17 точны при вычислении моментов до 11 и 17 порядков, соответственно, на всей сфере, а квадратуры ВОТ строились, исходя из требования, что моменты должны быть точны на главных полусферах. Тем не менее, таблица 3.1 показывает, что квадратуры Лебедева, вообще говоря, сравнимы по точности с квадратурами ВОТ при вычислении моментов на главных полусферах, чем квадратуры ВОТ. Более того, можно ожидать такого же поведения и при вычислении моментов, отличных от тех, что представлены в таблице 3.1.

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

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

Простейшая (вырожденная) «интерполяция» — «кусочно-постоянная». Квадратуры из [261] основаны фактически на такой «интерполяции». Однако кусочно-постоянное представление функций имеет низкий порядок точности, поэтому необходима интерполяция более высокого порядка.

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

Соболев построил в [38] простейшие квадратуры Ньютона-Котеса на сфере, инвариантные относительно группы 7. В этих квадратурах узловые точки образуют сетки двух типов, которые связаны с двумя типами триангуляции на сфере. Эти триангуляции устроены следующим образом. Рассмотрим октаэдр, вписанный в сферу, вершины которого лежат на координатных осях. Грань октаэдра рассматривается как базовый треугольник, см. рис. 3.2, 3.3. Каждая грань октаэдра (сторона базового треугольника) разбивается на N, N = 1, 2,..., равных интервалов, N — 1-й точкой. Вершины грани октаэдра (являющиеся также вершинами самого октаэдра) и точки, разбивающие грани октаэдра (его ребра, стороны базового треугольника), принадлежат треугольной сетке (на октаэдре), которая образована меньшими треугольниками. В 1-м типе триангуляции стороны меньших треугольников параллельны базовому, во 2-м типе — перпендикулярны, см. рис. 3.2, 3.3. Сетка на октаэдре образуется вершинами меньших треугольников. Примеры триангуляций обоих типов даны на рис. 3.2, 3.3. Триангуляция и соответствующая сетка на сфере получается проецированием таковых с октаэдра на сферу. Очевидно, что стороны сферических треугольников, образующих триангуляцию на сфере, принадлежат дугам большого круга. Заметим, что в триангуляции 2-го типа некоторые из меньших треугольников пересекают координатные плоскости и не принадлежат одному октанту. Отметим также, что разбиения в работе [261] совпадают с триангуляциями 1-го типа.

Квадратуры Гаусса и Маркова более точны, чем квадратуры Ньютона-Котеса,

Рисунок 3.2. Триангуляция 1-го типа, N = 2, 3.

Рисунок 3.3. Триангуляция 2-го типа, N = 2, 3. Штриховые линии — триангуляция, основанная на квадратуре Ь11, предложенной Лебедевым.

поскольку они используют для каждой свободной узловой точки все возможные параметры, на сфере это координаты узловой точки (два угла) и вес. В квадратурах Гаусса и Маркова узловые точки распределены на сфере более равномерно, чем в квадратурах Ньютона-Котеса. В некоторых квадратурах, предложенных Лебедевым, узловые точки образуют сетки, очень близкие с сетками в триангуляции 2-го типа. Например, узловые точки квадратуры Ь11 и триангуляции 2-го типа для N = 2 состоят из узловых точек 1-го, 2-го и 3-го классов и подмножества 5-го класса. Сами сетки отличаются лишь координатами узловых точек 5-го класса (единственные свободные точки в квадратуре). Для удобства вместо координат О этих точек мы указываем их проекции (х°п, у(, z0íi) на октаэдр, и представляем здесь одну из точек, остальные точки получаются при помощи преобразований группы 1 17). Для сетки 2-го типа (х°п,у°п, ) = (1/6,1/6, 2/3), а для квадратуры Лебедева Ь11 (х°п,у°п) = (1/5,1/5, 3/5). Триангуляция, основанная на точках квадратуры Лебедева Ь11 там, где она отличается от триангуляции 2-го типа, показана на рис. 3.3 штриховыми линиями. Таким образом, варьированием координат узловых точек 5-го класса можно строить триангуляции квази-второго типа. В связи с этим возникает вопрос: верно ли, что относительная разница между наибольшей и наименьшей площадью сферических треугольников достигает наименьшего значения

для триангуляции, основанной на квадратуре L11? Ответ отрицательный. Разница минимальна для триангуляции, в которой (х°п,у°0,z°n) « (0.2152449, 0.2152449, 0.5695102).

Следующая квадратура Лебедева, узловые точки которой близки к сетке 2-го типа для N = 3, это квадратура L17.

3.1.3 Кусочно-квазилинейная интерполяция, основанная на триангуляции 1-го типа

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

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

^ у 1 z

х = N У= N, Z= N,

где гх, гу, гz = -N,...,N, при этом имеет место одно из равенств

Nх| + |Ъу| + |iz| = N ± 1.

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

Шixiyiz = -Uixiyiz ,

где х = - х, у = - у, z = - z.

Узловые точки, соответствующие заданной триангуляции (для некоторого N), совпадают с вершинами сферических треугольников uixiyiz, то есть эти точки имеют координаты

(х\

Зу ,

z

где

ы + Цу | + hz | = n.

Всего на сфере имеется 4 N2 + 2 таких точек. Множество этих точек инвариантно относительно группы симметрии октаэдра, причем

= - П .

3x3y3z ]x]y]z

х

зХ + Зу + з 2

дальнейшем будет использоваться как последовательная нумерация индексом т так и тройная нумерация индексом ]х ]у]х.

Перед тем как объяснить, что такое кусочно-квазилинейная интерполяция на сфере, заметим, что вектору О = ( Дх,Ду ) на сфере соответствует вектор

(хЪ\ 1 /дл

уЪ

|Дх1 + |Ду I + I

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

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

м

1(г, п) « 1п(г)фп(П), (3.1)

|п|=1

где 1п(г) = 1(г, Оп), фп — элементарные кусочно-квазилинейные функции первого типа (конечные кусочно-квазилинейные элементы первого типа на сфере), определяемые равенством фп(От) = 8пт, 8пт символ Кронекера (то есть функция фп отлична от нуля только в тех треугольниках, для которых точка Оп является вершиной). В этом случае объемная плотность энергии излучения может быть аппроксимирована по формуле

» м

/ 1(г, п') « 4ж V тп1п(г),

Л2 ^

|п| = 1

(3.2)

где

1

т,п =

-1 [ фп(п')ап', (3.3)

а интеграл рассеяния может быть приближенно вычислен по формуле

1 Г м

— Ф(Пт, &)1(г, П')ап' ^^ттпШ, (3.4)

§2 I I 1

|п|=1

где

1

т

-1 [ Ф(Пт, &)фп(о')ая', (3.5)

} §2

|т|, |п| = 1,..., М. Если рассеяние изотропно, то есть Ф = 1, то ттп = тп.

Для вычисления интеграла рассеяния можно использовать также вместо форму-

лы (3.1) приближенное следующее представление подынтегральной функции

м

Ф(Пт, п'№, о') Ф(От, Оп) 1п(г)фп(П', (3.6)

|n| = 1

тогда интеграл рассеяния будет приближенно равен

м

1

4тг

I Ф(ПтП')1(г, П') dtf « V Ф(Пт, Qn)wJn(r).

Js |n|=1

т, аиn)Wn!n{r). (3.7)

Формула (3.7) имеет такой же вид, что и квадратурная формула для вычисления интеграла рассеяния в методе дискретных ординат.

Интегрирование по сфере сводится к интегрированию по сферическим треугольникам Шгхгугг. В параграфе 3.А.1 показано, как вычислять интегралы

[ /(п (3.8)

гхгугх

Интегрирование функций фп по сферическим треугольникам Uixiyiz описано в параграфе 3.A.2.

В таблице 3.2 представлены квадратурные веса wn для некоторых значений N. Относительная точность различных квадратурных схем при вычислении моментов

/ ПахПъуП1 dn

Jz>0

по полусфере z > 0 представлена в таблице 3.3. Квадратуры, основанные на весах wn, обозначаются PQLA^, N =3, 4,... (Piecewise QuasiLinear Angular, триангуляция 1-го типа). Эти квадратуры точны для моментов (0, 0,0) и (0, 0, 2).

Таблица 3.3 демонстрирует, что квадратуры PQLA^, вообще говоря, более точны, чем SN квадратуры, предложенные в работе [174], и менее точны, чем квадратуры, предложенные в работе [113]. Однако здесь следует заметить, что квадратуры PQLA разработаны для того, чтобы трактовать зеркальное рассеяние, а это не по силам другим квадратурам. Это преимущество квадратур PQLA перед другими отмечено в недавнем обзоре [88].

В следующем параграфе будет показано, что использование второго типа триангуляции позволяет существенно снизить ошибки в квадратурных формулах. Из таблицы 3.3 следует, что асимптотическая точность квадратурной схемы PQLA^, которая пропорциональна здесь 1/ М, М = 2N2 + 1, достигается для моментов, представленных в таблице, при N = 5.

Таблица 3.3 показывает также, что квадратуры Tn, которые можно интерпрети-

Таблица 3.2. Квадратурные веса 'шг,

N П = ]х]у 1 001 0.16667

2 002 0.02914

101 0.06876

3 003 0.01002

102 0.02675 111 0.03724

4 004 0.00492

103 0.01211

202 0.01696

112 0.01985

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

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

3.1.4 Кусочно-квазилинейная интерполяция, основанная на триангуляции 2-го типа

Кусочно-квазилинейная интерполяция, представленная в предыдущем параграфе, реализуема лишь для триангуляции 1-го типа. Для произвольной триангуляции, в частности, для триангуляции 2-го типа, кусочно-квазилинейная интерполяция может быть реализована следующим образом. Построим многогранник, вписанный в сферу, вершины которого совпадают с вершинами сферических треугольников. Пусть (хп, уп, ) — радиальная проекция точки О на многогранник. Непрерывную функцию /(О), О Е §2, будем называть кусочно-квазилинейной второго типа, если она линейна в координатах (хп, у^, г^) на каждой грани многогранника. Построение такой функции на сферическом треугольнике описано в параграфе З.А.З.

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

Таблица 3.3. Относительные ошибки при вычислении моментов /г>0 Щ по полусфере г > 0 для различных квадратурных схем.

Квад- Число Относительные ошибки (в %) для моментов (а, Ь, с)

ратура узлов

(0, 0, 1) (0, 0, 3) (0, 2, 1) (0, 0, 4) (0, 2, 2) (0, 0, 5) (0, 4, 1) (0, 2, 3) (2, 2, 1)

pqlaS,1' DCT020" DCT0206 L11 pqla22)1 38 48 48 50 50 -4.17 0 2.20 -1.23 -2.27 -0.14 0.28 -0.22 0 1.62 -8.20 -0.28 4.60 -2.48 -6.16 -0.68 0 0 0 3.05 1.03 0 0 0 -4.60 -0.73 -0.28 0.05 0 4.49 -14.32 -3.60 7.04 -2.48 -4.32 0.99 1.40 -0.80 0 -4.13 -8.25 6.46 8.07 -7.50 -15.78

pqla22)2 pqla22)2° pqla22)3 pqla41) T3 50 50 50 66 72 -0.66 -0.69 0 -2.48 0.10 -0.58 -0.56 -1.82 0.64 1.62 -0.76 -0.80 1.88 -5.60 -1.41 -1.12 -1.08 -3.55 1.30 3.24 1.70 1.63 5.30 -1.98 -4.85 -1.66 -1.57 -5.11 1.88 4.41 -1.04 -1.12 -1.04 -5.52 3.52 1.55 1.47 4.78 -1.81 -3.97 -4.48 -4.53 4.81 -13.57 -11.08

LSE8 LSO8 LSH8 EWE8 EWO8 88888 о о о о о 1.70 0 0 1.49 0 -0.14 0 0.22 -0.10 0 3.54 0 -0.22 3.08 0 0 0 0 0 -0.19 0 0 0 0 0.28 0.03 0.65 -0.19 0.02 0 5.52 -5.82 -1.24 4.68 -6.19 -0.47 -1.29 1.02 -0.35 0 5.61 20.07 0.39 5.17 18.56

LSN8 DCT111" DCT1116 pqla51) 80 80 80 80 102 2.07 0 0 1.30 -1.60 -0.95 0.82 0.02 -0.08 0.40 5.08 -0.82 -0.04 2.64 -3.60 -1.82 1.13 0 0 0.73 2.73 -1.70 0 0 -1.08 -2.86 1.07 0 0.02 0.98 6.36 2.01 -1.44 3.92 -3.76 2.89 0.33 -0.12 -0.24 -0.74 5.65 -0.82 3.87 4.50 -8.89

L17 pqla32)2° LSH10 EWE10 EWO10 110 110 120 120 120 -0.56 -0.42 0 1.02 0 0 -0.2 0.17 -0.05 0 -1.12 -0.64 -0.17 2.09 0 0 -0.4 0.10 0 0.01 0 0.58 -0.15 0 -0.01 0 -0.58 0 0.01 0 -1.44 -1.2 0.21 3.17 -0.32 0.01 0.57 0.52 -0.16 0 -2.49 -1.36 -2.72 3.36 0.97

T4 PQLA61) LSE12 S12 pqla71) 128 146 168 168 198 0.13 -1.12 0.97 1.44 -0.82 0.47 0.30 -0.04 -0.82 0.24 -0.21 -2.52 1.99 3.70 -1.88 0.84 0.55 0 -1.71 0.43 -1.26 -0.85 0 2.57 -0.63 1.16 0.74 0.05 -2.70 0.56 0.55 -2.56 3.05 4.31 -1.92 -0.91 -0.57 -0.14 2.93 -0.43 -1.07 -6.23 3.06 3.41 -4.64

T5 pqla81) ^16 T6 pqla91) 200 258 288 288 326 0.08 -0.64 1.10 0.05 -0.50 0.32 0.18 -0.65 0.22 0.14 -0.16 -1.44 2.84 -0.11 -1.16 0.63 0.33 -1.38 0.42 0.25 -0.94 -0.48 2.06 -0.64 -0.40 0.91 0.44 -2.16 0.61 0.35 0.48 -1.44 3.30 0.33 -1.12 -0.85 -0.34 2.37 -0.57 -0.28 -0.71 -3.59 2.43 -0.52 -2.85

T7 pqla10) T8 392 402 512 0.04 -0.41 0.03 0.16 0.12 0.12 -0.08 -0.92 -0.06 0.31 0.20 0.23 -0.46 -0.33 -0.35 0.45 0.29 0.34 0.24 -0.96 0.18 -0.42 -0.22 -0.32 -0.37 -2.32 -0.28

Точные значения ж ж/2 ж/4 2ж/5 2ж/15 ж/ 3 ж/8 ж/12 ж/ 2

Ьц и L17 квадратуры, Лебедев [23];

SN квадратуры, Lee [174];

LSOg квадратура, Lathrop and Carlson [172];

LSE^, LSO^, LSH^, EWEw и EWOw квадратуры, Fiveland [113]; LSN8 квадратура, Wakil and Sacadura [273]; TN квадратуры, Thurgood, et al. [261]; DCT квадратуры, Koch, et al. [160]

гралы вычисляются аналогично тому, как это выполнялось для триангуляции 1-го типа при помощи формул (3.1)—(3.7), но теперь фп — элементарные кусочно-квазилинейные функции второго типа (конечные кусочно-квазилинейные элементы второго типа на сфере). Интегрирование по сфере также сводится к интегрированию по сферическим треугольникам, аналогично тому как это описано в параграфе 3.A.2.

Приведем в качестве примера квадратуры, полученные при помощи кусочно-квазилинейной интерполяции 2-го типа. Эти квадратуры обозначаются PQLA^, N = 2, 3 (Piecewise QuasiLinear Angular, триангуляция 2-го типа), где значения I = 1, 2, 3 соответствуют различному расположению узловых точек.

(2)1

Квадратура PQLA2 соответствует оригинальной триангуляции 2-го типа из работы [38]:

w(1) ъ 0.00620996; w(2) ъ 0.02649607; w(3) ъ 0.02468077;

w(4) = 0; п ъ 0.2357023, wf] ъ 0.01894174; w(6) = 0.

42

В квадратуре PQLA22)2 узловые точки совпадают с узловыми точками в квадратуре Лебедева L11:

и 0.01053022; ^(2) и 0.02212533; ^(3) и 0.02249134; = 0; г 1 и 0.3015113, т{5) и 0.02041330; ^(6) = 0, Квадратура PQLA22)3 соответствует триангуляции, в которой относительная разница между наибольшей и наименьшей площадями сферических треугольников минимальна:

ш(1) и 0.01314671; 1^(2) ^ 0.01995346; 1^(3) ^ 0.02113951;

w(

(4) = 0; п ъ 0.3333216, w(5) ъ 0.02115905; w(6) = 0.

Таблица 3.4 сравнивает относительную точность квадратур PQLA^)г при вычислении моментов П^ПуП^ с№ по полусфере г > 0. Таблица демонстрирует, что расположение квадратурных узлов существенно влияет на точность квадратуры.

Таблица 3.4. Относительные ошибки при вычислении моментов /г>о ^х^у^ ^^ по полусфере г > 0.

Квадратура

pqla22)1 PQLa22)2 PQLA22)2a PQLa22)3

PQLA32)2a

Число узлов

Относительные ошибки (в %) для моментов (а, Ь, с)

50 50 50 50 110

Точные значения

(0,0, 1) (0, 0, 3) (0, 2, 1) (0, 0, 4) (0, 2, 2) (0, 0, 5) (0, 4, 1) (0, 2, 3) (2, 2, 1)

-2.27 1.62 -6.16 3.05 -4.60 4.49 -4.32 -4.13 -15.78

-0.66 -0.58 -0.76 -1.12 1.70 -1.66 -1.04 1.55 -4.48

-0.69 -0.56 -0.80 -1.08 1.63 -1.57 -1.12 1.47 -4.53

0 -1.82 1.88 -3.55 5.30 -5.11 -1.04 4.78 4.81

-0.42 -0.2 -0.64 -0.4 0.58 -0.58 -1.2 0.57 -1.36

ж ж/2 ж/4 2ж/5 2ж/15 ж/ 3 ж/8 ж/12 ж/ 2

С увеличением N расстояние между гранями многогранника и сферой уменьшается, и разница между плоскими и сферическими треугольниками становится прене-

брежимо мала. В этом случае интегрирование по сфере в интеграле (3.3) может быть аппроксимировано интегрированием по многограннику, и формула (3.3) принимает вид

1 (3.9)

12-к

где Sn — сумма площадей сферических треугольников, имеющих узловую точку Qп вершиной. Чтобы понять, какую ошибку вносит аппроксимация (3.9), мы строим квадратуры псевдо- pqla22)2 и псевдо- pqla32)2 , в которых веса вычисляются по формуле (3.9), эти квадратуры обозначаются PQLA22)2a and PQLA^2)2a, соответственно:

PQLA22)2a:

w(1) ъ 0.010628; w(2) ъ 0.022103; w(3) ъ 0.022464; wf] ъ 0.020410;

PQLA32)2a:

w(1) ъ 0.00376575; w(2) ъ 0.00996536; w(3) = 0; w(4) ъ 0.00962567;

w(5) ъ 0.00802249; w25) ъ 0.00973055; w35) ъ 0.01002474; w(6) = 0.

Сравнение квадратур pqla22)2 и pqla22)2" и таблица 3.4 показывают, что приближенная формула (3.9) аппроксимирует формулу (3.3) с высокой точностью уже при N = 2. Поэтому можно ожидать, что разница между квадратурой pqla32)2< , основанной на приближенных весах (3.9), и квадратурой pqla32)2 , основанной на точных весах (3.3), будет весьма мала. Кроме того, стоит отметить, что разница между весами квадратуры PQLA32)2a и квадратуры Лебедева L17 не превышает 1-2%.

Таблица 3.3 показывает, что квадратуры pqla22)2 и pqla32)2" весьма точны и сравнимы с лучшими квадратурами, имеющими близкое число узловых точек.

В заключение заметим, что плохая точность квадратур PQLA^ и PQLA22)1 указывает на весьма неоднородное расположение узловых точек в оригинальных квадратурах первого и второго типа в статье Соболева

3.1.5 Дискретизация граничного условия

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

1(г, П) = 1ехь(г, П)+ ра1 (г, П) + —[ \п • П'\1(г, П', п • П < 0, (3.10)

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

h = п - 2(П • п)п

направление луча, зеркально отраженного в направлении О, см. рис. 3.4.

п

5

Б

Сг ! с Рисунок 3.4. Зеркальное отражение.

После угловой дискретизации граничное условие (3.10) должно быть записано в узловых точках От, п • От < 0 (узловые точки нумеруются последовательно так, что |га| = 1,... , М и О-т = — От, то есть число узловых точек равно 2М). С использованием представления интенсивности в виде (3.1) эта процедура осуществляется элементарно. Подставляя представление (3.1) в граничное условие (3.10), полагая О = От, и учитывая, что 1(г, От) = 1т(г), мы получаем граничное условие

м м

1т(г) = 1ехь(г, Пт) + Рз ^ 1п(г)фпфт) + рЛ ^ 1п(г)ап(г), П • Пт < 0, (3.11)

Н = 1 |п| = 1

где

1

ап(г) = —1 I п • П'фп(П') ¿.О'

К Зп-П' >0

К .)п-О'>0 и

Пт = — 2(От • п)п.

Коэффициенты ап всегда удовлетворяют тождеству

м

]>>п = 1, (3.12)

|п|=1

которое выражает закон сохранения энергии. Заметим, что коэффициенты ат в граничном условии (3.11) играют ту же роль, что и сумма У]3=1 щ1т,мт при вычислении потока излучения на границе по квадратурной формуле в оригинальном методе дискретных ординат (пг и Iт, г = 1, 2, 3, — координаты нормали п и направляющие косинусы направления 17т, соответственно).

Граничное условие (3.11) это дискретный аналог граничного условия (3.10). Кажется, что первое содержит формальное противоречие. Действительно, суммирование в (3.11) выполняется по всем направлениям Оп, в то время как в условии (3.10) учитываются только направления, удовлетворяющие условию п • Оп > 0 (падающее на границу излучение). В действительности (большая) часть направлений, для которых п • Оп < 0, не дает никакого вклада в правую часть условия (3.11) поскольку соот-

ветствующие коэффициенты фп(Ит) и ап(г) равны нулю. Тем не менее некоторые из этих коэффициентов, соответствующих таким направлениям, отличны от нуля: это направления, имеющие по крайней мере одно смежное направление Оп/, удовлетворяющее условию п Оп! > 0), то есть близкие к касательным направлениям. В результате правая часть условия (3.11) может зависеть от интенсивностей 1п(г), для которых п • Оп < 0, что с физической точки зрения бессмысленно. Чтобы преодолеть эту бессмысленность, можно исключить соответствующие узловые точки из суммирования, однако такой прием может привести к значительным ошибкам. Лучший способ состоит в экстраполяции интенсивностей, в направлениях, близких к касательным. Простейшая экстраполяция может быть реализована следующим образом. Пусть Оп — направление, для которого п • Оп < 0 и коэффициенты фп(£2т) или ап(г) отличны от нуля. Тогда интенсивность 1п(г) в правой части условия (3.11) может быть экстраполирована усреднением интенсивностей 1п/ (г) по всем направлениям Оп/, смежным с Оп и удовлетворяющим условию п • Оп! > 0. Более общая линейная экстраполяция также может быть легко реализована.

В заключение заметим, что для вычисления интеграла в условии (3.10) можно использовать вместо представления (3.1) приближенное представление подынтегральной функции

м

п • п' 1(г, п') « ^ п • Оп1п(г)фп(0'), (3.13)

|п| = 1

тогда вместо последней суммы в условии (3.11) будет

м

- р^п •ПпЪ(г)Ьп(г), (3.14)

|п|=1

где

Мг) = -/ фп(&)ЛП'. (3.15)

3.2 Уравнение переноса излучения и граничные условия в случае осевой симметрии

Этот раздел основан на результатах статей [228, 230]. 3.2.1 Уравнение переноса

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

стационарным уравнением переноса излучения [29, 32, 143, 196]

d¿ + üd¿ + Qdi 'дх у ду z dz

ü •VI + pl = пх— + Üv— + nz— +PI = asSI + к!ъ, (3.16)

где I = I(r, ü) (или в другой записи I(r, в, )) — интенсивность излучения,

r = (г cos р, г sin р, z)

— координата точки в пространстве (здесь сразу удобно перейти к цилиндрическим координатам),

ü = (sin в cos рп, sin в sin рп, cos в)

— направление излучения, к и as — коэффициенты поглощения и рассеяния, соответственно,

3 = к + as

коэффициент ослабления,

SI = SI(r) = [ I(r,ü)dü,

4TT J S2

— интеграл рассеяния, d ü = sin (9 d в dрп (для простоты рассеяние предполагается изотропным),

т - Т( ) п2(тТА

4 = 4(r) =-

ж

— интенсивность излучения черного тела, Т = Т(r) — температура, п — показатель преломления среды, т — постоянная Стефана-Больцмана.

Если задача для уравнения переноса осесимметрична (осью симметрии является ось z), интенсивность излучения I зависит только от переменных r,z,6 и |р — рп|. Поэтому достаточно искать значения интенсивности только в точках

(r, в, рп = 0) = (г cos р, г sin р, z, в, рп = 0), р Е [0, ж]

(углы р Е (ж, 2ж) исключены из соображений симметрии), то есть искать решение уравнения переноса (3.16) в половине {у ^ 0} трехмерной цилиндрической области и только для направлений

ü = (sin в, 0, cos 0), (3.17)

см. рис. 3.5.

В результате получаем, что уравнение переноса в случае осевой симметрии при-

У

f Ü

1 Divm \ x

Рисунок 3.5. При наличии осевой симметрии уравнение переноса достаточно рассматривать в половине {у > 0} трехмерной осесимметричной области D и лишь для направлений П = (sin в, 0, cos 0).

нимает вид

di di

ü • V I + 31 = sin б1— + cos б1— + 31 = asSI + к1ь, dx dz

(3.18)

3.19)

где I = I (г, в), интеграл рассеяния вычисляется по формуле

1 f^ f^

SI = SI(r, z) = — I(r, 0)sin0d edp,

J о Jo

r = (r cos p, r sin p, z). Решение уравнения переноса (3.18) ищется в области

(г, 0) eD{y^o} х [0,Н где B{y>0} — половина осесимметричной области D, для которой у > 0, см. рис. 3.5.

3.2.2 Другой способ вывода уравнения переноса

Этот способ вывода уравнения переноса (3.18) поможет вывести для него граничные условия.

В случае осевой симметрии стационарное уравнение переноса излучения в цилиндрических координатах имеет вид [29]

sin

cos p-

d (г I) d (sin p I)

d d p где I = I(r, z, в, p) — интенсивность излучения,

P = pn -'p>

d

+ cos 9——+31 = a„SI + к Ih, dz

(3.20)

(3.21)

pr — азимутальный угол в цилиндрических координатах, определяющий (вместе с г и ) координату

г = (г cos pr, г sin pr, z)

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

V

ляющий (вместе с полярным углом ) направление

ü = (sin в cos рп, sin в sin рп, cos в)

излучения в точке r, к — коэффициент поглощения,

т - т( ) п2тТ4

h = h(r) =-

ж

— интенсивность излучения черного тела, Т = Т(r) — температура, п — показатель преломления среды, ( — постоянная Стефана-Больцмана.

Благодаря соотношению (3.21) уравнение переноса (3.20) может быть интерпретироваться двумя способами. С одной стороны, полагая рг = 0 (и, следовательно, р = pf¡), получаем, что решение уравнения переноса (3.20) ищется в точках

r = ( , 0, )

для всех направлений

ü = (sin б1 cos р, sin б1 sin р, cos в), р = рп.

С другой стороны, полагая рп = 0 (и, следовательно, р = — рг), получаем, что решение уравнения переноса (3.20) ищется в точках

r = (г cos р, г sin р, z), р Е [0,ж], (3.22)

(то есть в половине цилиндрической области, в которой углы р Е (ж, 2ж) исключены из соображений симметрии), но только для направлений ü, удовлетворяющих условию рп = 0, то есть (3.17). Отметим, что в статьях [6, 7, 86] такой переход был фактически проделан, но ни уравнение (3.18), ни граничные условия для него получены не были.

Для применения метода характеристик второй способ интерпретации удобнее и естественнее, поскольку в этом случае мы имеем дело с обычным трехмерным физическим пространством, в котором характеристики являются прямыми линиями. Кроме того, координаты ( , р, ) могут быть интерпретированы как обычные цилиндрические координаты. Преобразовав уравнение переноса (3.20) к декартовым координатам (х, у, z), где х = г cos р, у = г sin р, получим уравнение (3.18), где I = I(r, в), координата точки r определяется формулой (3.22), а полярный угол в задает направление ü (3.17) (здесь интерпретация независимых переменных в целом (в совокупности) отличается от их интерпретации в исходном уравнении переноса (3.20)).

Заметим, что в то время как характеристики уравнения (3.20) криволинейны: г sin p = r0 sin p0, z — z0 = ctgO (r cos p — r0 cos p0), характеристики уравнения (3.18) — прямые линии:

У = Уo, z — zo = ctg0(;r — жо).

3.2.3 Условие на непрозрачной диффузной границе

Выведем граничные условия для уравнения переноса (3.18). Рассмотрим непрозрачную диффузно отражающую границу. Граничное условие для уравнения (3.20) имеет вид