Книга
← К списку книгЗагрузка разделов…
активнаАктивные записи участвуют в каталоге и синках. Неактивная книга скрыта, но не удалена из базы.
Книга
Предисловие автора
Одно из направлений науки XX века получило название тепло-массо-перенос. Развитие реактивной авиации, космических ракет, ядерной энергетики, атомного оружия, в частности, связано с успехами в описании тепло-массо-переноса. Высадиться на Луну, спроектировать атомную подводную лодку или водородную бомбу невозможно без предварительного численного моделирования. В число рассчитываемых физических процессов входят энергоперенос излучением и газодинамика. Рекордные по сложности вычисления тепло-массо-переноса выполняются на ЭВМ.
Численным моделированием физических процессов занимаются коллективы физиков и математиков. У нас в стране принято говорить о математической физике и вычислительной математике. Американцы говорят о вычислительной физике. Направление данной работы, с точки зрения автора, можно характеризовать как математическую и алгоритмическую физику.
В 1958-2002 гг. автор работал в Российском Федеральном Ядерном Центре (г. Саров) во Всероссийском научно-исследовательском институте экспериментальной физики (ВНИИЭФ). Посчастливилось участвовать в освоении почти всех отечественных ЭВМ, в освоении и развитии ряда методик и программ счета тепло-массо-переноса. В 1977-1997 гг. автор возглавлял лабораторию расчетов переноса излучения в математическом отделении ВНИИЭФ. В узкий круг интересов автора входило моделирование нестационарного переноса излучения с использованием метода коэффициентов видимости. Метод является рекордсменом по точности и экономичности численного решения двумерных нестационарных краевых задач энергопереноса излучением внутри ограниченной полости сложной формы с частично отражающими стенками.
Во ВНИИЭФ сложился авторский коллектив физиков и математиков, совместно работавших над уточнением описания кинетических эффектов в рамках расчетов энергопереноса излучением в диффузионном приближении. В нее входили, например, чл.-корр. АН СССР Ю. Н. Бабаев, Н.А. Дмитриев, В. Г. Заграфов, В. Е. Трощиев, С. В. Баженов, Б.П. Тихомиров, П.А. Перепелкин, П.И. Певная. Успеху работы способствова-
ли тесные контакты с Институтом прикладной математики АН СССР (г. Москва, акад. И.М. Гельфанд, Н. Н. Ченцов и др.) и с Всероссийским научно-исследовательским институтом прикладной физики (г. Снежинск, акад. Е.Н. Аврорин, М.Ю. Коз- манов и др.). Во ВНИИЭФ получены и основные излагаемые здесь результаты.
В ограниченной области пространства решается краевая задача нестационарного энергопереноса излучением в рамках корпускулярного подхода. Волновые и релятивистские эффекты переноса излучения не описываются. Описываются, в частности, большинство явлений переноса излучения в природе, наблюдаемых нами как говорят, невооруженным глазом.
На корпускулярном подходе основаны три эквивалентные модели описания переноса излучения: фотометрическая модель, интегро-дифференциальное уравнение и интегральное уравнение Пайерлса. Известные численные дискретные модели корпускулярного переноса лучистой энергии основаны на двух аналитических формах описания переноса: интегро-дифференциального уравнении и интегральном уравнении Пайерлса. На каждом шаге по времени многократно итерациями совместно решается полная система аппроксимирующих алгебраических уравнений для всех счетных источников излучения и всех счетных направлений переноса. Итерации направлены на согласование различных видов описания энергии для достижения с приемлемой точностью общего баланса энергии и получения возможности счета нового шага по времени. Используются потоковые методы описания переноса. Применяются интерполяции, требующие от решения повышенной гладкости по сдвиговым переменным.
В 1969 г. по инициативе автора была начата разработка нового подхода к численному решению двумерных и трехмерных задач на основе фотометрической модели. Фотометрическая модель была первой теоретической моделью, построенной на обобщении опытных данных на основе корпускулярного подхода. В число опытных данных входило предположение о некогерентности источников излучения.
Следуя фотометрической модели, предполагалось получить численную модель переноса с независимым рассмотрением точечных источников излучения в силу предположения их некогерентности. Стремились избежать трудности решения системы аппроксимирующих уравнений.
Планируемые цели достигнуты на минимальном теоретическом и алгоритмическом уровне. Излагаемые результаты работы получены в два этапа.
На первом этапе работы на основе нового подхода удалось с 1972 г. решать нестационарные по излучению задачи внутри неподвижных областей сложной формы с пренебрежимо малым взаимодействием излучения и вещества. Граничные поверхности поглощают и частично отражают излучение. Учитывается запаздывание и затенение. Счет идет по явной балансной схеме с контролируемой точностью. При счете двумерных задач существует возможность счета с заданной высокой точностью, обеспеченной полуаналитическим расчетом коэффициентов видимости.
В творческую группу вместе с автором вошли Р. Ф. Маши- нин, Э.Н. Тихомирова. Новая программа стала первой программой для решения трехмерных задач нестационарного энергопереноса. С шестью точными значащими цифрами в разнообразных вариантах постановки граничных условий решается, в частности, трехмерная стационарная тестовая задача В. Г. Рогачева о переносе излучения внутри сферической полости, имеющая аналитическое решение.
Результаты работы были опубликованы спустя 12 лет после начала счета задач. Было решено до тысячи трехмерных задач. Автором всех этих расчетов была Э. Н. Тихомирова.
Второй этап работ выполнен в 1999-2001 гг. в соавторстве с Е.А. Карповцевым, Е. В. Морозовой, И. А. Нарожной, В. А. Но- вичихиным, Э.Н. Тихомировой и Е.И. Рябовым. Удалось перейти к описанию переноса с учетом взаимодействия излучения и вещества. Использовались модификации написанных ранее модулей — в соавторстве с Э. Н. Тихомировой, Р. Ф. Машининым и А. И. Кирилловым. С шестью точными значащими цифрами была решена тестовая задача — известная нестационарная задача Р. Д. Рихтмайера и К. Мортона о выходе импульса излучения из центра шара с учетом изотропного рассеяния.
После этого предельным переходом было получено аппроксимируемое уравнение — новая аналитическая форма интегрального уравнения переноса излучения.
С учетом запаздывания и затенения описывается распределение импульса точечного источника:
— угловое распределение, распределение вдоль расходящегося пучка лучей,
распределение выведенной веществом части энергии внутри объема расчетной области,
распределение проникающей части излучения на граничной поверхности.
Новый дискретный подход на основе метода коэффициентов видимости обеспечивает явное описание распределения импуль-
са излучения точечного источника в соответствии с фотометрической моделью с учетом взаимодействия излучения и вещества. Расчет некоторого числа шагов по времени выполняется в два этапа:
расчет коэффициентов распределения,
расчет переноса энергии на основе вычисленных коэффициентов.
Аналитически коэффициенты описываются двукратными определенными интегралами (интегралами по поверхности граничных ячеек). Вычисление значения подынтегральной функции включает пространственные геометрические расчеты пересечения фиксированного луча (из точки источника в точку- цель на границе) с поверхностью ячеек внутренней сети. Геометрия расчетной области определяется, вообще говоря, в ходе газодинамического расчета, совместного с расчетом переноса излучения.
Сумма коэффициентов распределения для фиксированного источника должна равняться единице. Отличие суммы вычисленных коэффициентов от единицы — средство контроля локальной точности счета. Нормирование значений коэффициентов делением на сумму обеспечивает абсолютную балансность счета.
Расчет переноса — циклы умножения на коэффициенты — описание распределения энергии очередного источника по всем сеточным целям.
Новый характеристический подход к решению нестационарных задач кажется слишком сложным. Большинство математиков и физиков привыкло к описанию переноса на шаге по времени лишь в соседние с источником ячейки. В новом методе отсутствуют вычисления, не имеющие ясного физического смысла, есть полнота описания лучевого переноса на всю глубину расчетной области — на каждом шаге по времени.
По определению Р. Беллмана и Р. Калабы новая схема счета при решении двумерных задач близка к квазианалитичности: вычисления сведены к процедурам, выполняемым с заданной точностью. Это вычисление двукратных определенных интегралов и итерации Ньютона для решения алгебраических уравнений четвертой степени. В каждой внутренней ячейке с помощью итераций Ньютона решается уравнение локального баланса поглощенной, внутренней и излученной энергии.
Новая схема счета расщепляется по источникам и целям, она удобна для параллельных вычислений.
Автор с благодарностью вспоминает о многолетней коллективной работе. Каждый из ее участников внес свой творческий
вклад, без которого не было бы общего успеха. Особенно автор благодарен Е.А. Карповцеву, ставшему организатором последнего цикла работ. Автор выражает благодарность сотрудникам РФЯЦ-ВНИИЭФ, РФЯЦ-ВНИИТФ и ИПМ РАН за обсуждение результатов работы и их публикацию в сборнике статей ВАНТ. Автор признателен В.Ф. Юдинцеву, В. А. Жмайло, А. Н. Субботину, А. В. Никифоровой, Н.П. Плетеневой и И.М. Анисиной — за критические замечания и помощь в работе, А. В. Забродину, И. Ф. Подливаеву, Б. П. Тихомирову и В. В. Гореву за поддержку в работе, И.Д. Софронову за прекрасные условия работы в руководимом им творческом коллективе математиков ВНИИЭФ. Работа велась в тесном контакте с заказчиками расчетов, физиками-теоретиками ВНИИЭФ. Открытый обмен информацией сдерживался лишь секретностью расчетов. Все вместе способствовало развитию методик и программ счета.
Введение
Настоящая публикация посвящена изложению теоретических и алгоритмических особенностей нового лучевого подхода к решению двумерных и трехмерных краевых задач нестационарного энергопереноса излучением [1-35] на основе фотометрической модели и алгебраических законов теплопереноса в условиях локального равновесия излучения и вещества. Учитывается взаимодействие излучения и вещества с сохранением класса гладкости решения, с квазианалитическим [9] вычислением сеточных коэффициентов распределения. Учитывается запаздывание и затенение. Изложение ограничивается простейшими аналитическими моделями описания углового распределения излучения точечных источников, приближением серой материи, случаем изотропного рассеяния.
Мы не ставим целью полноту обзора литературы по рассматриваемой тематике, она обширна. Разработка физических и численных моделей переноса лучистой энергии в корпускулярном приближении является одной из актуальных тем развития экспериментальной и математической физики. Терминология и обозначения наиболее близки к [7, 8, 10, 11]. Часть обозначений будем определять при первом обращении к ним, поясним наиболее часто встречающиеся в тексте термины.
Со времен И. Ньютона и Г. Лейбница описание процессов, связанных с излучением, базируется на двух подходах, корпускулярном и волновом. Первая стадия корпускулярного описания переноса излучения завершилась в XVIII веке построением опытно-теоретической фотометрической модели [1, 22, 25, 32].
Большую роль при формировании этой модели сыграл опыт наблюдения переноса зрительной информации через малое отверстие в стене темной комнаты (камеры-обскуры). Бесконечное число лучей, идущих от разных точек, не смешиваясь между собой, проектируют изображение на экран, переворачивая его.
Этот опыт свидетельствовал о перспективности использования понятия точечного источника излучения, о большой угловой делимости излучения точечного источника, о корпускулах — частицах, переносящих изучение вдоль пучка расходящихся лучей. Опыт свидетельствовал о возможности на этой, казалось бы
бесконечно малой, основе получить вторичные точечные источники, и, более того, о возможности повторить перенос еще раз, чтобы пронести информацию от вторичных источников к глазу наблюдателя.
Другие опыты привели к необходимости рассмотрения другой, волновой, модели переноса излучения.
Основным для фотометрической корпускулярной модели стало предположение о некогерентности точечных источников излучения, дающее возможность их независимого рассмотрения. Энергия, излученная разными источниками в одном направлении, складывается. Лучи от разных источников, имеющие разные направления, не взаимодействуют друг с другом. Фиксированный луч взаимодействует с веществом среды переноса. Совместное влияние некогерентных источников приводит к изменению состояния вещества и к изменению состояния всех точечных источников — суммарному результату независимого переноса лучистой энергии от изолированных источников. Теоретическое описание взаимодействия излучения и вещества было разработано позже завершения разработки фотометрической модели.
Точечный источник, внутренний или граничный, излучает энергию внутрь полного телесного угла или его части. Эта энергия распределяется по расходящимся лучам. Авторами фотометрической модели были разработаны ставшие классическими аналитические варианты описания углового распределения и переноса излучения вдоль луча. Фотометрическую модель переноса излучения замкнули закон Ламберта, дающий пример аналитического описания углового распределения энергии для точечного источника на плоской поверхности изотропно излучающего вещества и закон Бугера-Ламберта аналитического описания ослабления интенсивности излучения при прохождении луча сквозь полупрозрачный слой вещества.
Современная стадия корпускулярного описания переноса излучения началась в конце XIX — начале XX веков. Сначала были разработаны уравнения Больцмана, описывающие перенос широкого класса взаимодействующих частиц. Для частиц со свойствами корпускул излучения математически следовало интегро-дифференциальное уравнение переноса излучения. Математическим же путем на основе интегрирования этого уравнения в первой половине XX века было получено интегральное уравнение Пайерлса переноса излучения.
В главах 1, 2 мы изложим основные сведения об известных уравнениях переноса, интегро-дифференциальном и интегральном уравнении Пайерлса, которые будут использованы в гл. 3 при
выводе интегрального уравнения в новой форме, форме описания распределения лучистой энергии точечного источника. Изложенный материал глав 1-3, в основном, следует тексту публикаций [34, 35]. В гл. 3 включены материалы, полученные автором в процессе работы с рукописью данной публикации. Глава 4 содержит описание алгоритмов первой схемы численного решения двумерных нестационарных задач на основе нового уравнения с помощью метода коэффициентов видимости. Глава 5 содержит описание алгоритмов вычисления коэффициентов. Глава 6 посвящена краткому обзору решения трехмерных задач энергопереноса излучением внутри неподвижных полостей сложной формы [18]. Глава 7 излагает первые результаты решения двумерных задач с учетом взаимодействия излучения и вещества [33]. Изложенный материал, в основном, следует тексту публикаций [18, 28, 33].
Частицы излучения, корпускулы, при описании переноса имеют разные названия, например, фотоны, кванты. К настоящему времени два последних термина слились, фотон — название элементарной частицы — кванта поля электромагнитного излучения. В нашем изложении частицы излучения выступают как доли энергии, которые, в свою очередь, можно делить на части и складывать. Предполагается, что каждая описываемая доля энергии содержит большое число квантов энергии. Доли энергии описываются кусочно-непрерывными функциями, вопрос о кратности этих долей числу квантов носит условный характер.
На основе фотометрической модели мы будем строить схему аналитического и численного приближенного решения нестационарных краевых задач переноса излучения. Рассматривается фиксированный точечный источник — точка с прилегающей окрестностью границы или объема.
За некоторый промежуток времени, dt или At, следующий за фиксированным моментом времени t, точечный источник излучает некоторую порцию энергии. Будем называть эти кратковременные повторяющиеся порции энергии — импульсами излучения точечных источников. Аналитическое описание распределения импульсов фиксированного точечного источника является новой формой интегрального уравнения переноса излучения.
Известное интегральное уравнение переноса излучения — уравнение Пайерлса — описывает перенос излучения всех точечных источников в фиксированную точку-цель вдоль пучка сходящихся лучей.
Излучение фиксированного точечного источника, распределенное по расходящимся лучам, пронизывает вещество расчетной области, вступая с ним во взаимодействие. Проникающая
часть излучения достигает границы области и распределяется по ней. Выведенная веществом часть излучения распределяется внутри объема расчетной области. Доли распределенного импульса, достигшие цели — элемента объема или границы, начинают складываться с такими же долями распределенной энергии других источников. Рассчитывается плотность одностороннего потока энергии, проникающей на границу, и плотность выведенной энергии внутри объема. На замыкающей стадии вычислений, относящихся к одному промежутку времени излучения, dt или At, — локально для каждого элемента граничной поверхности или объема, рассчитывается новый импульс энергии излучения на новом промежутке времени.
Используются законы Стефана-Больцмана и Кирхгофа, закон Ламберта углового распределения граничного излучения, закон Бугера-Ламберта ослабления излучения при переносе вдоль луча, принцип Шустера разделения величины выведенной веществом энергии на слагаемые величины поглощенной и рассеянной энергии. Локально для каждого элемента объема вещества итерациями Ньютона решается нелинейное уравнение баланса энергии: часть поглощенной энергии идет на повышение внутренней энергии и температуры вещества, часть — на повышение плотности собственного излучения вещества.
Расчетная область V объема вещества ограничена кусочногладкой граничной поверхностью S. Внутри расчетной области в виде заданных функций определены параметры вещества и начальные значения счетных величин. Гладкость функций определяется возможностью их приближения кусочно-постоянными сеточными функциями. На граничной поверхности заданы граничные условия. Внутри области и на ее границе могут быть определены участки и функции заданного ввода лучистой энергии. Расчет энергопереноса ведется в некотором диапазоне времени t. Численные расчеты выполняются с постоянным шагом At по времени, At = т.
Рассматриваемый метод объединяет сеточное определение параметров вещества, точечных источников, дискретных пучков расходящихся лучей, и геометрический расчет распределения энергии излучения точечных источников вдоль лучей.
Фиксированные точки будем обозначать большими буквами Р, Q, R, сообщая этим их принадлежность объему или границе и роль в описании переноса — точка-источник, точка-цель.
Геометрическим местом фиксированного точечного источника может служить либо точка Р границы, Р Е S, тогда говорят о граничном точечном источнике Р, или, еще короче, о граничной
точке-источнике Р. Либо им может служить точка Р объема, Р EV, тогда говорят о внутреннем точечном источнике Р или о внутренней точке-источнике Р. Численные расчеты выполняются для точек-источников Р, являющихся счетными центрами ячеек, граничных или внутренних. Как правило, точка Q — граничная точка-цель, конец вектора PQ описания переноса. Точка R — внутренняя точка элемента объема, пересекаемого лучом.
С помощью рассуждений, носящих предварительный характер, попробуем пояснить суть подхода к построению излагаемой численной схемы описания переноса.
Пусть объем расчетной области разбит на ячейки внутренней сети, и граничная поверхность разбита на ячейки граничной сети. Пусть известны координаты точечных источников Р внутри каждой ячейки внутренней и граничной сетей. И пусть известны правила дополнительного разбиения граничных ячеек на составляющие элементы и правила вычисления координат точек-целей Q внутри каждого элемента дополнительного разбиения. Тогда определены лучи PQ из фиксированного точечного источника в точки-цели на границе. Дополнительное разбиение граничных ячеек на составляющие элементы используется лишь на стадии расчета коэффициентов распределения — для приближенного вычисления определенных интегралов по телесному углу. На стадии численного описания переноса происходит распределение импульса энергии для каждого сеточного источника в отдельности с помощью умножения на коэффициенты распределения, в сумме составляющие единицу.
В каждый рассчитываемый момент времени t (на шаге по времени) нестационарного расчета для каждого точечного источника Р вычисляется энергия Е^РД), излучаемая им за предстоящий промежуток времени т (величина шага по времени). Рассматривается дискретный пучок расходящихся лучей PQ с центром в точке-источнике Р и точками Q цели на элементах дополнительного разбиения граничных ячеек, Q Е S. Рассматриваются отнесенные к лучам PQ доли ^^р^ телесного угла, опирающегося на элементы дополнительного разбиения граничных ячеек. Для дискретного пучка лучей PQ и долей ^^р^ телесного угла рассчитывается угловое распределение импульса энергии Е(Р, t\
По пути от источника до цели каждый луч пересекает ячейки объема области V, рассчитывается распределение вдоль луча PQ выведенной веществом части энергии. Требуется обозначе-
ние для промежуточныхточек R, Re PQ, R Е V. Пусть распределение вдоль луча PQ выведенной веществом энергии вычисляется для пересекаемых лучом внутренних ячеек. Тогда за множество точек R можно принять множество счетных центров ячеек внутренней сети. До конечных граничных точек Q доходит проникающая часть излучения.
Запишем сказанное в виде формулы для случая конечного числа лучей PQ в пучке и конечного числа промежуточных точек R на луче —
Е(РИ) =^E(P,t)pPQ + £ E(P,t)pPRQ. (1) Q R,Q
Размерность левой части (1): [Е(РД)] = Вт. Первое слагаемое правой части (1) описывает распределение на границе проникающей части излучения. Вдоль каждого дискретного луча PQ до граничной точки Q дойдет некоторая доля ppq первоначальной энергии Е(Р, t\
Второе слагаемое правой части (1) описывает распределение выведенного веществом (поглощенного и рассеянного) излучения. На каждом дискретном луче PQ мы фиксируем доли Pprq первоначальной энергии Е^РД\ выведенные веществом в окрестности промежуточных точек R объема V на этом луче.
Продолжая далее рассуждения, носящие предварительный характер, будем избегать введения временных обозначений для промежуточных величин.
Пусть множество граничных точек Q теперь является множеством счетных центров не элементов разбиения граничных ячеек, а самих ячеек. Так можно будет считать после объединения суммированием порций энергии, вынесенных разными лучами на все элементы дополнительного разбиения каждой фиксированной ячейки граничной поверхности. При этом объединятся соответствующие слагаемые в первой сумме правой части (1). Теперь доля ppq — доля импульса энергии Е(РД) источника Р, достигшая при переносе граничной ячейки с центром Q.
Пусть множество промежуточных точек R — множество счетных центров ячеек внутренней сети. Может случиться, что одна и та же промежуточная точка R встречается во второй сумме правой части (1) несколько раз — для разных лучей, и соответствующие доли Pprq выведенной энергии относятся к разным лучам PQ, но к одной ячейке внутренней сети с центром R. Обозначим ррр сумму долей Pprq для фиксированной промежуточной точки R, но для разных точек Q цели, т. е. для разных
лучей переноса (доли pprq определены до объединения долей энергии, вынесенной разными лучами в одну и ту же граничную ячейку), PPR = ^PPRQ.
Q
Теперь доля ррр — доля импульса, исходящего из источника Р, достигшая при переносе внутренней ячейки с центром R. Тогда (1) перепишется в виде
Е(р,#) = 52 Е(р, ^PPQ + 52 Е(р, t^PPR’ (2) Q R
'pPPQ A-'pPFR= 1. (3) Q R
Обнаруживается характерная особенность описания распределения энергии точечного источника на основе фотометрической модели переноса. Имеется два варианта формул распределения энергии точечного источника.
В одном варианте, (2), описывается распределение конкретного импульса. В левой части равенства — величина излученного импульса. Правая часть формулы содержит два слагаемых. Описываются доли распределения по ячейкам границы проникающего излучения и доли распределения внутри объема выведенного веществом излучения. Все члены уравнения содержат множителем величину импульса Е^РД).
После сокращения всех членов (2) на величину импульса Е(РД) получаем второй вариант формулы описания переноса энергии в виде равенства (3) для величин долей распределения: в сумме они равны единице.
Формула (3) явится основой для вычисления коэффициентов распределения. Формула (2) служит для расчета переноса лучистой энергии.
Подчеркнем, что множество лучей первоначально не связано с множеством граничных ячеек. Чем больше мы будем рассматривать лучей, первоначально ведущих в элементы дополнительного разбиения граничных ячеек, тем лучше мы будем описывать заполнение лучами объема. Тем точнее мы будем описывать перенос с учетом взаимодействия излучения фиксированного источника с веществом ячеек объема и вынос проникающего излучения в граничные ячейки. С помощью дополнительного разбиения можно менять число лучей, добиваясь пересчетом нужной точности описания углового распределения энергии для фиксированного источника.
Мы описали распределение импульса, излученной точечным источником, по ячейкам граничной и внутренней сети. Часть лучистой энергии точечного источника Р, выведенную веществом ячеек внутренней сети, с учетом запаздывания опишем формулой Ep^RR + pPR/с) = pprE^PR). Проникающую часть лучистой энергии, распределенную по ячейкам граничной сети, с учетом запаздывания опишем формулой Ep(Q,t + ppq/c) = ppqE(PR). Введено обозначение для счетных частиц (2) распределения излучения фиксированного точечного источника — для частиц, которыми ячейка-источник с центром в точке Р обменивается с ячейками-целями с центрами Q, R.
С учетом запаздывания фиксировано время прихода частиц, Ppq, Ppr — средние расстояния до целей, с — скорость света,
Е(Рр) = ^Ep^Q^t + ppq/c) + ^Ep^RR Р Ppr/с).(4)
QR
Введем обозначения для целых, rpq, прр, и дробных, apQ, арр, частей деления расстояний на путь света за шаг по времени,
Ppq = (rpq + «fq)ct, ppp = (npR + аРКут.(5)
Импульс E(P,t) и доли распределения этой энергией, Ep(Q,tf), Ep(R,tff), связаны с диапазоном шага по времени. Если бы все расстояния Ppq, Ppr были кратны пути ст света за шаг по времени, мы бы имели вместо (4):
Е(Р, t) = у EP(Q, t + npQ-т) + У Ep(R, t + пррт) (6) QR
Через несколько, npQ, шагов по времени пришедшая в граничную ячейку Sq доля энергии фиксированного источника Р станет долей Ep(Q,tf) одностороннего потока энергии. Вся плотность потока в момент времени tr будет вычисляться суммированием доли энергии от фиксированного источника Р с другими одновременно приходящими долями энергии — от каких-то других источников (не обязательно излучившими свой импульс энергии в момент времени t — из-за разных расстояний).
Пришедшая во внутреннюю ячейку Vr доля энергии фиксированного источника Р через несколько шагов по времени станет долей Ер(РДгг) выведенной веществом этой ячейки энергии. Вся выведенная энергия ячейки в момент времени tn будет вычисляться суммированием с другими одновременно приходящими долями выведенной энергии — от каких-то других источников.
Анализ на основе линейной модели излучения энергии в промежутке времени, dt или At, приводит к следующему результату.
В общем случае, когда расстояния ppg, ррр не обязательно кратны пути ст света за шаг по времени, обменные порции энергии Ep(Q,trY Ep(R,trr) приходят в ячейки-цели на двух последовательных шагах по времени своими долями
(1 - apQ^Ep^Q, t + tipqtY aPQEp{Q,t + ripQT + т),
(1 - app)Ep(R,t + прртУ aPpEP(R,t + пРрт + т\
Сумма долей ppq, ррр равна единице, (3). Ясен физический смысл этих дискретных долей-вероятностей, которые мы будем называть коэффициентами распределения энергии фиксированного источника вместе с числами tipq, aRQ, и прр, арр. Доли Ppq — вероятности для излучения фиксированного источника Р попасть в телесные углы лучей PQ и дойти без поглощения и рассеяния до ячеек-целей Q граничной сети. Доли ррр — вероятности поглотиться или рассеяться веществом ячеек внутренней сети. Числа tirq, aRQ, и ripp, арр, (5), позволяют согласовать перенос с учетом дискретных шагов т по времени.
При непрерывном описании суммы должны перейти в интегралы, а дискретные коэффициенты распределения энергии Ppq, Ррр — в функции распределения энергии Krq, Крр.
Мы прервем рассуждения, носящие предварительный характер и поясняющие интегральное описание переноса энергии точечного источника излучения на основе фотометрической модели.
Два ряда обозначений, например dt или At = г, вводятся в связи с рассмотрением в дальнейшем двух вариантов описания переноса излучения: аналитического и численного дискретного. Особенностью излагаемого здесь подхода к описанию переноса излучения является близость этих двух вариантов описания. Большинство счетных формул численного метода являются дискретными вариантами аналитических формул.
Для указания окрестности фиксированных значений независимых переменных при численном описании нам нужны конечные дискретные величины, например At = г; при аналитическом описании нужно рассмотрение бесконечно малых положительных величин, например dt, в процессе стягивания дискретных окрестностей к фиксированным значениям независимых переменных, например t.
Обозначим окрестности семи независимых переменных:
спектральный интервал, div или Az/, частоты z/;
окрестность объема, dVp или Vp, вокруг внутренней точки Р расчетной области — по всем трем координатам;
окрестность площади, dSp или Sp, вокруг точки Р границы — по двум координатам;
элемент телесного угла, da)p^ или Aw-^, окрестность вектора PQ — конус с вершиной Р;
диапазон времени, dt или At = г, для момента времени t.
Будем рассматривать точки Р граничной поверхности лишь на гладких участках поверхности S, предполагая, что для точки Р и для ее открытой окрестности dsp определена касательная плоскость. Будем предполагать, что точкам линий излома граничной поверхности будет соответствовать множество меры ноль — по площади прилегающей окрестности. Излучение всей граничной поверхности будет описываться точечными источниками на гладких участках поверхности. При дискретном рассмотрении линии излома граничной поверхности будем включать в множество линий разбиения граничной поверхности на счетные ячейки. При этом мы, в частности, учтем разрыв решения в точках излома границы, [10].
Переходим к описанию точечных источников. Нестационарный точечный источник — это точка Р со своей излучающей окрестностью и момент времени t и его окрестность, шаг по времени. Точка Р определяется своими координатами (множество меры ноль). Окрестность этой точки — окрестность объема, dVp или Vp, излучающего вещества для внутренней точки или окрестность площади, dSp или Sp, поверхности излучающего вещества для точки границы. В фиксированный момент времени в течение некоторого промежутка времени, dt или г, для окрестности объема или границы, точечный источник Р излучает лучистую энергию. Испускается импульс (порция) энергии, dE(P,t) или E(P,ty размерности \dE(P,t), E(P,t)] = Вт:
для внутреннего точечного источника
dE(P, t) = ит,8(Р, t)dVp dt или Е(Р, t) = иТ)8(Р, t)Vpr, (7)
для граничного точечного источника
dE(P,t) = J+(P,t)dSPdt или Е(РД) = Jp(t)SPr. (8)
Ниже мы поясним физический смысл функций иГ)8(РДУ J+(P,t) и других рассчитываемых функций. Из (7) и (8) ясно, что функции иГ)81 J+ описывают распределение излучательной способности вещества. В точке Р в момент времени t определяется энергия, излучаемая веществом единицы объема в единицу времени. Размерность [nrs(F, t)] = Вт/(м3с), U4P^ = Вт/(м2с).’
Сначала об основных функциях J^P^t), J+(P,t), описывающих излучение на границе S расчетной области V, и о граничных условиях нестационарной краевой задачи, с помощью которых замыкаются вычисления граничных величин — для одного шага по времени.
При решении краевых задач переноса в граничных точках- целях рассчитывается величина J~(Q,tf) плотности одностороннего потока проникающего излучения — от разных точечных источников, [2-5, 7]. Верхний индекс (—) соответствует направлению, обратному к направлению вектора внутренней нормали пр. Внутреннее направление нормали пр к граничной поверхности S в точке Р границы определяется по отношению к расчетной области V.
При описании точечных источников для граничной поверхности рассчитывается величина Р^РД) плотности обратного одностороннего потока лучистой энергии, высвечиваемой с границы внутрь расчетной области, [2-5, 7]. Верхний индекс (+) соответствует направлению вектора внутренней нормали пр в точке Р к граничной поверхности S. Размерность плотности односторонних потоков [J-, J+] = Вт/(м2с).
По определению обе функции J-, J+ могут принимать лишь неотрицательные значения.
Предполагается, что плотности односторонних потоков J+, J- связаны между собой граничным условием
P(F,f) = /p(J"(F,f)),(9)
где функция /р задана для каждой точки Р границы и позволяет вычислить J+ на основе J”. Функция /р позволяет вычислить плотность потока, высвечиваемого внутрь расчетной области на основе плотности потока, вычисленного в ходе описания переноса. Мы получаем возможность рассчитать величину импульса энергии для граничных точечных источников на новом шаге по времени, (8). Вычисления граничных величин для одного шага по времени замкнулись.
Приведем два примера граничных условий (9).
На некоторых участках границы может быть задана функция альбедо а(РД) — отношение односторонних потоков,
J+(P,t) = a(P,t)J~ (P,t)-(Ю)
Внутрь расчетной области высвечивается заданная доля энергии, проникающей на границу при переносе вдоль лучей от точечных источников.
С помощью граничного условия (10) не всегда удобно описывать ввод энергии с границы. Удобнее использовать определение плотности q(P, t) потока энергии через границу — в окружающее пространство. На некоторых участках границы может быть задана функция q(P, t) — разность односторонних потоков,
q^P,t)=J-(P,t)-J^P,t),
J~(P,t) = J+(P,t) + q(P,t)- (И) Односторонний поток энергии, проникающей на границу, частично высвечивается обратно внутрь расчетной области, частично уходит в окружающее пространство. По определению поток q(P,t) в окружающее пространство определен как положительная функция, если односторонний поток направлен против вектора внутренней нормали, т. е. если J > JP Если задан ввод энергии в область расчета в виде потока q^Pp), то функция q^Pp) должна быть отрицательной.
Описание источников излучения, (7), (8), использует функции J^(P,ty иг^РЛ\ рассчитанные на предыдущем шаге. Описание переноса завершается вычислением функций J-(P,t), Ua,s(P,t\
Перейдем к описанию основных функций для точек объема вещества, заданных и рассчитываемых. По определению все функции расчета, которые будут введены ниже, могут принимать лишь неотрицательные значения.
Следуя [7, гл. 2], рассмотрим сначала функцию распределения энергии излучения /(z/,F, Q,t), функцию семи независимых переменных: частота излучения iv, три координаты точки Р расчетной области V, два угла наклона к осевому направлению для фиксирования вектора FQ переноса излучения, наконец, время t. Определением функции f^v.P.Pip) служит выражение:
f(iv,P,Pl,t)divdVPdu^ (И)
— для числа квантов спектрального интервала div относительно частоты iv — внутри окрестности объема dVp вокруг внутренней точки Р, внутри элемента телесного угла du— около вектора переноса FQ, в момент времени t.
Тогда формула
UUU= | (12)
(4тг) 0
определяет плотность энергии излучения всего спектрального диапазона частот — энергию единицы объема — вокруг внутренней точки Р расчетной области в момент времени t.
Здесь hv — квант энергии частоты р, интегрирование выполнено для всех частот и для всех направлений переноса.
Из сравнения (11) и (12) может показаться необычным переход к единичному объему — замена единицей бесконечно-малой величины dVp. Здесь и в дальнейшем нужно придерживаться правила математического анализа: действия с конкретными бесконечно малыми величинами правильны, если правильны эти же действия с любыми определяющими их конечными малыми величинами в конкретном предельном процессе. Если, например, в расчете известна энергия /\Ер для счетной ячейки с объемом АУр, то плотность энергии U(P) вычисляется по формуле U(P) = ^ЕР!АУр. При переходе к бесконечно-малым величинам энергии 6ЕР и объема размещения этой энергии 5Vp — в процессе стягивания ячеек к точке Р — остается понятным и правильным смысл выражения U(P) = 8Ep(8Vp. Более того именно это выражение и определяет правильное значение плотности энергии в точке Р. Если мы будем пытаться искать плотность U(P) в точке Р как отношение каких-либо конечных величин, вычисленных в том же предельном процессе, мы, вообще говоря, будем получать разные значения для плотности при разных размерах ячеек АУр. Ведь ячейки АУр включают в себя и другие точки, не только точку Р, а в других точках может быть и другая плотность энергии. Бесконечно-малые величины 8ЕР и 6Vp являются бесконечно-малыми величинами одного порядка малости, а в таком случае использовать их отношение не только можно, как мы видим, но и нужно. Все непрерывные вычисляемые функции — физические величины — в таких аналитических описаниях, как решение краевой задачи, определяются на множествах меры ноль: в точках, для лучей — с помощью бесконечно-малых объемов, площадей, телесных углов, диапазонов времени. Численно значения этих величин устанавливаются в дискретных расчетах и уточняются в упомянутом выше предельном процессе повторения расчетов для измельченных вариантов сети ячеек. Вообще говоря, разные значения отношения АЕ^/АУр должны быть приближенными значениями для плотности U^P) и сходиться к точному значению в предельном процессе. После сказанного, возвращаясь к формулам (11) и (12), подчеркнем: формула (12) может быть получена заменой единицей бесконечно-малой величины dVp, ведь слева конечная величина U^P) относится к точке, но определяет энергию для единицы объема.
Читая текст, все время нужно прослеживать физический смысл формул и проверять размерность величин.
Функция П(РД) одномоментно — в момент времени t — описывает, (12), в единице объема около внутренней точки Р расчетной области совокупную энергию всех квантов излучения — частиц энергии, несущихся со скоростью света в этот момент времени внутри этой единицы объема. Энергия описывается для единицы объема, но определяется, исходя из состояния излучения внутри бесконечно-малой окрестности dVp точки Р.
Обозначим Vp какую-нибудь окрестность точки Р с объемом, равным единице.
Зная величину U(P,tY ничего нельзя сказать о состоянии излучения в окрестности dVq некоторой другой точки Q, даже если Q е КД dVQ с КД
Мы должны научиться обращаться с функциями распределения массы или энергии типа функции U(P,t) или объемной плотности вещества р^РД), описывающими состояние вещества в единице объема, но отражающими свойства вещества лишь в бесконечно малой окрестности точки Р и момента времени t. Опираться надо на физический смысл математических действий, полезен контроль размерностей слагаемых и суммы, множителей и произведения.
При описании переноса излучения с учетом взаимодействия излучения и вещества, мы будем придерживаться правила раздельного описания излучения с указанием места и причины его рождения. Внутренние источники излучения рождают кванты в процессе рассеяния и в процессе собственного излучения вещества в условиях локального термодинамического равновесия. При взаимодействии излучения и вещества часть излучения, выведенная веществом за счет поглощения и рассеяния, гибнет, — рассеянное излучение меняет направление. Вещество вновь рождает излучение изотропно — в рассматриваемом нами приближении — высвечивает рассеянную им энергию и энергию собственного излучения — на основе законов Кирхгофа и Стефана-Больцмана.
От одной независимой переменной v мы избавимся рассмотрением переноса интегрально по всем частотам — в приближении серой материи, как в черно-белом кино. Еще говорят, — в одногрупповом приближении. И новую форму интегрального уравнения, и соответствующую численную схему можно изложить и в многогрупповом спектральном приближении. Один из трехмерных расчетов переноса излучения внутри полости сложной формы был выполнен для тридцати спектральных групп.
Кроме функции плотности энергии излучения и^РД\ мы будем использовать неотрицательные функции состояния вещества и фракций плотности энергии излучения, связанных с описанием взаимодействия излучения и вещества,
И^РДУ И^РЛУ Xa_s(P,tY Ue(P,ty,
Ur(P,ty U.8YP,tY «r,s(F, ^) = ur(FF) + us(F, t),(13)
ua(P, t\ ua^s{P, t) = иа^Р, t) + u8(P, t\
Размерность этих функций: \иа, и8, иаД = 1/м, \U, Ue\ = = Вт/(м3), [ur, us, ua] = Вт/(м3с).
Нижние индексы поясняют физический смысл различных функций, описывающих взаимодействие излучения и вещества:
т — для описания испускательной способности вещества за счет собственного излучения в условиях локального термодинамического равновесия,
s — для описания испускательной способности вещества за счет рассеяния излучения,
а — для описания способности вещества к поглощению излучения,
е — для описания внутренней энергии вещества.
Заданные функции состояния вещества xaix8lxai8 — коэффициенты взаимодействия излучения и вещества: коэффициенты, соответственно, поглощения, рассеяния и полного взаимодействия; ха^8 = ха + х8. Коэффициент поглощения ха — заданная функция двух независимых параметров состояния вещества, температуры Т^РЛ^ и плотности р^РД); х^РД) = f(T,p\ Коэффициент рассеяния х8 — заданная функция плотности, xs(P,f) = f(p\ Кванты поглощаются и излучаются при переходах связанных электронов из одного энергетического состояния атомной системы вещества в другое, [7, гл. 2], при переходах связанных электронов с одной стационарной орбиты на другую. При этом происходит изменение температуры вещества, электрон переходит на другую стационарную орбиту или покидает атом, происходит ионизация вещества. Кванты рассеиваются, в нашей модели, на свободных электронах частично или полностью ионизованного вещества. Квант меняет направление. Квант с прежним направлением как бы погибает. Рождение кванта с новым направлением приписывается источнику излучения вещества, расположенного в окрестности гибели прежнего кванта. При условии изотропности рассеивания можно рассматривать суммарно энергию квантов, рассеянных веществом, и энергию квантов, излученных веществом, в предположении локального
равновесия. Собственное излучение вещества, в нашей модели, изотропно.
Функции Ue, uTl us, иа, иГ)8, uai8 — основные вычисляемые функции, описывающие излучение и взаимодействие излучения и вещества во внутренних точках Р расчетной области в момент времени t.
Например, иТ)8(РД) характеризует излучательную способность точечного источника Р в момент времени t: описывается суммарная энергия, излученная единицей объема вещества в единицу времени — за счет собственного излучения (г) и за счет рассеяния (s);
функция иаЦРД) описывает суммарную энергию, выведенную единицей объема вещества в единицу времени — за счет поглощения излучения (а) и за счет рассеяния излучения (s) — выведенную из всех пересекающих лучей описания переноса в соответствии с коэффициентами поглощения и рассеяния и законом Бугера-Ламберта.
Левая часть интегро-дифференциального уравнения описывает перенос излучения с помощью дифференциального оператора и функции интенсивности, 1ЦР, O,t),
ЦР, n, t) = JIvdv = c^ M<V, Р, П, t) dv,0 0
| ЦР, Q, t) du = c | | hvf dvdu = cU(P, t). (4тг) (4тг) О
Функция интенсивности описывает угловое распределение плотности потока энергии точечного источника и распределение энергии вдоль расходящегося пучка лучей. В фиксированный момент времени описывается распределение энергии излучения фиксированной частоты — в точках пространства и по направлениям. От одной независимой переменной v мы избавимся рассмотрением переноса интегрально по всем частотам — в приближении серой материи. Размерность интенсивности ЦР, О,£) в приближении серой материи: [1(F, fi,t)] = Вт/(м2сср), это плотность направленного потока энергии — величина энергии, переносимой из точки Р в направлении вектора РО в единичном телесном угле (ср — стерадиан) вокруг этого вектора, через единицу площади, перпендикулярной вектору О, в единицу времени.
На основании (14) обсудим физический смысл функций U, cU, I. Умножение на скорость с света плотности лучистой энергии U приводит к описанию плотности потока энергии точечного внутреннего источника Р внутрь полного телесного угла.
Для изотропных источников излучения, когда функция распределения / и интенсивности Д не зависят от направления переноса О, плотность излучения равна
Uv(P,t) = 4Trhvf =^IV, (15)
после интегрирования (15) по частоте
П(РД) = ^cU(P,t) = I(P,t). (16)
Одинакова размерность \^cU(р, t) Л^ЦРЛ^ = Вт/(м2с). При изотропности излучения источника Р, (16) определяет направленный поток I в любом направлении, он одинаков, не зависит от направления.
В плотность U в (16) входят только кванты, родившиеся в точке Р изотропного источника — за счет рассеяния и за счет собственного излучения — и не входят кванты, родившиеся в других точках и перенесенные в окрестность точки Р.
На основе опытных данных, полученных французским ученым П. Бугером в 1729 г., немецкий ученый И. Г. Ламберт в 1760 г. разработал первую теоретическую модель одномерного переноса излучения вдоль луча, ведущего от точечного источника,
dl = -xldg. (17)
За счет взаимодействия с веществом излучение убывает. Дифференциал интенсивности, описывающий вывод энергии излучения веществом, пропорционален самой интенсивности излучения (ее называли яркостью) и дифференциалу пути излучения, dg = cdt, от источника. Коэффициент пропорциональности х является характеристикой вещества и называется коэффициентом взаимодействия излучения и вещества.
В нашем изложении мы будем выделять два физических процесса взаимодействия: поглощение и рассеяние, х = xQvS = = ха + xs. С учетом (16) из (17) следует: выведенная по всем направлениям, ведущим от источника, энергия пропорциональна плотности излучения,
(18)(18)^a,s — ^a,s^U, Us — XscU Ua — XacU.
(18)
(18)
Пропорциональность в (18) следует из (17) и, в общем случае для выведенной по всем направлениям энергии — с учетом (14).
Выведенная энергия uai8 в нашей методике вычисляется непосредственно — из интегрального уравнения переноса излучения в новой форме. Методики, использующие аппроксимацию интегро-дифференциального уравнения и интегрального уравнения Пайерлса, определяют выведенную энергию на основе ее пропорциональности (18) плотности излучения. Плотность излучения U в остальных методиках является основной расчетной величиной, определяющей точность расчета. В нашем случае цикл вычислений вообще обходится без вычисления плотности излучения U. Плотность излучения может быть вычислена дополнительно в момент выдачи результатов счета, — например, с помощью (18).
Пусть Су — заданный параметр вещества, удельный коэффициент теплоемкости вещества, и вычислены значения двух независимых параметров состояния вещества, температуры Т^РД) и плотности р(РД\ Тогда определим функцию Пе как линейную функцию температуры, Ue = СурТ. В условиях локального равновесия собственное излучение вещества пропорционально четвертой степени температуры (законы Кирхгофа и Стефана- Больцмана),
ит = ^acUT, Ur = аТ4, (19)
а — коэффициент Стефана.
Функция UT в теории переноса излучения называется плотностью излучения, соответствующего точечному источнику Р абсолютно черного излучателя.
Мы кончили предварительное обсуждение основных функций описания переноса излучения. Нужно замкнуть расчет в точках объема.
Во внутренних точках R описание переноса начинается с вычисления энергии для точечных источников — на основе функции ur,s(R,t\ а приводит описание переноса к определению величины энергии, выведенной веществом окружающего точку R объема (IVr, и к определению функции НаДРД). С помощью принципа Шустера по сумме иа^8 = иа-V и8 с помощью отношений коэффициентов взаимодействия находят слагаемые иа, и8: Ua = ua,sxalxa,s, us = uai8xslxa,s. Как и в других методиках, см., например [8], при решении краевых задач для замыкания вычислений во внутренних точках локально для каждой точки решается нелинейное алгебраическое уравнение баланса энергии для вещества: часть поглощенной энергии идет на повышение
внутренней энергии и температуры вещества, часть — на повышение плотности собственного излучения вещества:
dUe
7,— ^а ^г-
at
При численной реализации дифференциальное уравнение заменяется разностным уравнением — алгебраическим уравнением четвертой степени относительно температуры вещества:
Ue(t + г) - Ue(t)CvpT(t + r)-ue(t)4
тт
(20)
Здесь Ue^ и иа — известные величины. Нелинейное алгебраическое уравнение четвертой степени решается итерациями Ньютона относительно неизвестной температуры Т.
Вычисленные локально для каждой внутренней ячейки функции uT)S = ur + us определяют внутренние точечные источники на следующем шаге по времени.
Мы замкнули расчет в точках объема.
Функция uT)S определяет правую часть интегро-дифференциального уравнения.
Продолжим рассуждения, носящие предварительный характер и поясняющие интегральное описание переноса энергии точечного источника излучения на основе фотометрической модели.
В дискретном случае мы будем описывать распределение энергии излучения точечного источника в пространстве и во времени, т. е. будем превращать сумму в слагаемые с помощью умножения суммы на числа — доли единицы — дискретные вероятности.
На аналитическом уровне описания распределения в виде интегрального уравнения в левой части равенства мы предпочитаем начинать с распределяемой, уже вычисленной, величины энергии одного точечного источника. А в правой части интегрального уравнения — с суммы двух интегралов, описывающих распределение фиксированного импульса одного источника по граничной поверхности и внутри объема расчетной области. Это первая форма нового интегрального уравнения — описания распределения энергии точечного источника излучения.
Вторая форма нового интегрального уравнения переноса излучения точечного источника проще первой формы, использующей сумму двух интегралов, как и интегральное уравнение Пайерлса. Выводу уравнения во второй форме посвящены главы 1 и 3. Интеграл по телесному углу описывает угловое распре-
деление энергии точечного источника. Внутри этого интеграла с помощью одномерного определенного интеграла описывается перенос излучения вдоль луча.
Вернемся к обсуждению на предварительном уровне особенностей описания распределения излучения точечного источника — к обсуждению первой формы соответствующего интегрального уравнения.
Интеграл по граничной поверхности описывает вклад фиксированного точечного источника в односторонний поток энергии, перенесенной на границу сквозь вещество области с учетом запаздывания и затенения — проникающая часть энергии излучения. Интеграл по объему описывает вклад фиксированного точечного источника в плотность выведенной веществом энергии излучения. Указывается время запаздывания для каждой доли распределенной энергии. Интегрирование производится по видимой из точки-источника части граничной поверхности и объема расчетной области. Кроме понимания физической картины распределения на этой стадии описания фотометрической модели в виде интегрального уравнения ничего не требуется.
С помощью (4)-(6) и (8) запишем сказанное в виде формул для граничного точечного источника
dE(P,t) = J+(P,t)dSPdt,
dE(P,t) = j dEp(Q,t + ppQ/c)dSQ +
S(F)
+ j dEp(R, t-V Ppr/c) dVp, (21) y(F)
(P, t) dSp dt = | Jp(Q, t + npQr) dSq dt +
+ | Ua^p^RH-VnpRT^dVpdt.
Уравнения математической физики, уравнения переноса излучения в том числе, описывают физические процессы на языке математического анализа — с помощью рассмотрения бесконечно-малых величин. Описываются физические связи между рассчитываемыми и заданными функциями состояния излучения и вещества, относящимися к единице площади или объема, — в определенных классах гладкости этих функций, например, интенсивность излучения, плотность вещества. И в то же время это
описание ведется с помощью рассмотрения бесконечно-малых величин одного порядка малости, бесконечно-малых окрестностей площади и объема — окрестности некоторых фиксированных точек — и бесконечно-малых окрестностей телесного угла — окрестности некоторых лучей описания переноса.
При аналитическом описании распределения излучения точечного источника мы неизбежно должны заниматься распределением бесконечно-малых величин, т. е. должны прибегать к рассмотрению бесконечно-малых величин разного порядка малости. Выручает простота опытно-теоретической фотометрической модели и следующая из нее наглядность геометрической оптики. Нужно настроиться на неизбежность рассмотрения непривычных аналитических конструкций, иных, чем рассматривались в известных описаниях переноса излучения, — непривычных, но целиком относящихся к математическому анализу.
В (21) слева — бесконечно-малая величина dE энергии точечного источника излучения, выраженная через конечную величину плотности одностороннего потока. Справа — описание распределения этого импульса бесконечно-малыми долями dEp — в точках граничной поверхности и объема. Ясно, что эти доли распределения импульса энергии — бесконечно-малые величины более высокого порядка, чем сам импульс. Но для рассчитываемых долей распределения тоже возможно их описание через плотности потока и плотность энергии,
dEp(Q, i) = Jp (Q, i) dSq dt, dEp(R, i) = ua,sp(R, t) dVp dt.
Функции одностороннего потока Jp, проникающего на границу излучения, и плотности ua,sp, выведенной веществом энергии — для одного фиксированного источника Р — описывают эти величины для единицы площади и объема — так же, как функции J”, uaiS для всех источников вместе. А по отношению к величинам одностороннего потока J“(Q,t) и плотности u^^R,^ — величины Jp, uasp — бесконечно-малые величины. Величины Jp, ^a,sP описывают вклад одного точечного источника — в односторонний поток J-(Q,t) и плотность uafS(R,t), соответственно. Этот вклад точечного источника обеспечивается с помощью рассмотрения одного луча и его окрестности телесного угла. Казалось бы, можно отразить это с помощью традиционного обозначения типа dJ, dua,s. Появление новых букв увеличивает трудность чтения, использование выражения dJ~dSqdt наряду с J~ (Q,tf) dSq dt ничего не добавляет информационно, кроме сообщения о порядке малости, наконец, величина dJ относится
к единице площади и времени. Поэтому мы ввели другое обозначение малости величины — для величин вклада в плотность потока и в плотность энергии — с помощью добавления нижнего индекса, Jp, ua,sp, — сообщения о точке источника. Полная величина рассчитываемых величин потока J- или плотности uaiS получается суммированием вкладов от всех источников (конец гл. 3).
Описание распределения именно импульса энергии позволяет избежать ошибок, связанных с использованием величин разной размерности. Как правило, потом требуется сокращение левой и правой частей полученного интегрального уравнения на бесконечно малую величину dt времени высвечивания импульса энергии, получается уравнение распределения мощности точечного источника,
J+(P,t)dSp= | Jp(Q, t + npqr) dSg +
+ | и^р^КДРтгррт^Ур.
V(P)
Мощность излучения точечного источника внутрь расчетной области в результате переноса вдоль пучка расходящихся лучей распределяется по граничной поверхности — в виде мощности одностороннего потока проникающего излучения — и распределяется внутри объема вещества — в виде мощности поглощения и рассеяния. Размерность членов уравнения, после сокращения на dt, [j+dSp, JpdSq, ur,sdVp, ua,spdVp] = Вт/с.
Метод контроля правильности действий с бесконечно малыми величинами известен — проверка этих действий с допредельными дискретными величинами, опять же, опираясь на физический смысл. Формулам предшествовали их дискретные аналоги. В дальнейшем изложении, наоборот, все интегральные действия предшествуют тем же действиям в дискретной форме. Обычно вычислительные математики так и поступают — строят численную модель с помощью физико-математической модели — дифференциальных или интегральных уравнений, описывающих физические процессы на основе эксперимента.
Изложенное содержит противоречие. Аналитические средства описания нагляднее, обозримее, чем дискретные, но требуют необходимой подготовки от читателя, известной степени свободы в специфических, как правило, аналитических преобразованиях с величинами, операторами и функционалами. Нередко требуется рассмотрение новых непривычных сначала аналитических
средств, без которых невозможно изложить физическую модель. Так, например, была введена в рассмотрение 5-функция Дирака.
мы ввелимы ввелиНаш материал проще, аналитические средства и физический смысл знакомы студентам, физикам и математикам, технических вузов и университетов. Но это первое аналитическое и алгоритмическое описание фотометрической модели лучистого переноса излучения точечного источника. Автор не считает, что это первое изложение не нуждается в совершенствовании. Первое изложение как раз направлено на то, чтобы привлечь внимание к новому подходу. В данном случае, — к новой возможности сделать шаг к квазианалитичности вычислений при решении нестационарных краевых задач переноса излучения, прежде всего, — в полупрозрачных областях со сложной геометрией границ с учетом движения вещества, взаимодействия излучения и вещества, с учетом запаздывания и затенения.
мы ввели
мы ввели
точечноготочечногоВ дискретном варианте описания распределения обозначение для долей ppq, ppp импульса Е(РД) источника излучения — см. (2), (3),
точечного
точечного
(22)(22)е(р, г) = £ ppqPp, tHE pprE^P t)QR
(22)
(22)
E ppQ + E ppR = 1 •QR
При непрерывном описании суммы должны перейти в интегралы, а дискретные вероятности ppq,ppr — в функции плотности распределения энергии KpQexp (—tpq), КррехрРтрр), аналитическое описание которых нам предстоит в дальнейшем,
dE(P, f) = dE(P, f) j Kpq exp (-tpq) dSq +
+ dE(P,t) | KpRexp(-TpR)dVR, (23)
dE(P, t) = dE(P, t) j KPq exp PtPq^ dSQ +
+ dE(P, t)E j KPR exP Ptpr) EVr.
Vr VR(P)
Подынтегральные функции (23) описывают распределение величины импульса. Мы разбили граничную поверхность S и объем V на ячейки Sq, Vr, как в (22). При этом области видимости 8^Р\ V(P) разбились на области видимости для ячеек, Sq^P), Vr(P\ Мы получили интегральное описание для коэффициентов распределения при численном решении краевой задачи,
dE(P, t) = dE(P, t) £ ppQ + dE(P, t} £ ppR,QR
Ppq= j Кр9ехр(-трУ^89,(24)
Sq(P)
Ppr= I Kppexp (-tpr) dVp.
Ясен физический смысл функций распределения Krq, Krr и известна основа для аналитического описания этих функций и вычисления после этого коэффициентов распределения (24). Функция KPQdSQ для излучения граничного источника должна являться функцией распределения трех независимых событий. Нужно описать распределение энергии по направлениям PQ. Нужно описать распределение энергии по элементам dwpg телесного угла, составляющим окрестность лучей PQ, ведущих на границу. Нужно описать ослабление для пучка dupQ излучения на пути до границы. Примеры соответствующих функций распределения известны, их аналитический вид используется в интегральном уравнении Пайерлса. Это закон Ламберта углового распределения граничного излучения, формула (28), см. ниже, — для определения доли телесного угла и закон Бугера-Ламберта ослабления излучения — при описании переноса вдоль луча,
KpQdSQ = cos (пР, PQWpa = dSQ,' ^PpQ
Рр = cos (пр, PQ},rq = cos (ng, QF),
I KpQ exp (-tpq} dSQ = SqWI KpQ exp (-tpq} dSQ = SqWj KPq exp (-tPq) dSQ = | exp (-tpq} cos (nP,PQ}dWp^, SqWSQ{P)
I KpQ exp (-tpq} dSQ = SqW
I KpQ exp (-tpq} dSQ = SqW
^P-tpq^^MSq. (25)
<PPQ
SqW
2 Ю. А. Дементьев
Множитель exp (—тру) описывает ослабление потока лучистой энергии вследствие взаимодействия излучения с веществом на пути PQ — закон Бугера-Ламберта; описывается доля энергии проникающего излучения по отношению к первоначально излученной энергии; tpq — оптическая длина пути PQ,
Qtpq = | ^a,s dg.р
Множитель др — косинус угла наклона вектора PQ к вектору пр внутренней нормали к граничной поверхности в точке Р источника; вероятность излучиться в направлении PQ для источника излучающей граничной поверхности пропорциональна др — закон Ламберта. Мы будем предполагать, что граничная поверхность в окрестности dSp точечного источника Р излучает по закону Ламберта и для нее в точке Р определена касательная плоскость.
Выражение
d.sQ
Ppq
описывает проекцию площади dSq на сферу единичного радиуса с центром в точке Р источника излучения.
Проекционная геометрия тесно связана с геометрической оптикой, описанием пучка расходящихся из точки лучей. Рассматривается сфера единичного радиуса с центром в точке Р и телесные углы — конические поверхности с той же вершиной. Телесный угол измеряется площадью пересечения конической поверхности со сферой единичного радиуса — объемный аналог линейного угла на плоскости, измеряемого длиной дуги окружности единичного радиуса.
Ниже будет показано, что интегрирование выражения dujpQ
по полусфере телесного угла приводит к единице, можно это выражение принять за функцию распределения энергии граничного точечного источника на стадии описания углового распределения энергии. Ниже будет показано, что интегрирование выражения
по сфере телесного угла приводит к единице, можно это выражение принять за функцию распределения энергии внутреннего
точечного источника на стадии описания углового распределения энергии.
Пусть площадь dSg — площадь окрестности граничной точки-цели Q, часть площади ячейки Sq граничной поверхности. Мы примем телесный угол dcopQ,
dujpq = -^- dSg,
Ppq
за телесный угол — окрестность луча PQ при описании распределения энергии точечного источника Р. Это проекция площади dSq на сферу единичного радиуса.
После определения аналитического описания функции Kpq можно приближенно вычислить интегралы — коэффициенты распределения ppq — и рассчитать распределение проникающей части излучения — для граничного и внутреннего точечного источника.
С функцией Крр сложнее. В связи с приведенным выше рассуждением для дискретного случая можно ожидать, что функция распределения энергии для части излучения, выведенной веществом, определена для фиксированного луча. Вновь опорой будет служить закон Бугера-Ламберта ослабления излучения — при описании переноса вдоль луча.
Мы завершили рассмотрения, носящие предварительный характер, выражения (1)-(25), и поясняющие интегральное описание переноса энергии точечного источника излучения на основе фотометрической модели.
Основные аналитические формулы в нашем изложении, содержащие интегралы и описывающие распределение импульса энергии точечного источника с помощью пучка расходящихся лучей, это формулы прямого действия в наиболее точном выражении — по отношению к дискретным аналогам этих формул. Определенные интегралы пополнили список аналитических функций. Они отнесены к квазианалитическим средствам вычисления, см. [9].
При закрепленной сети граничных и внутренних ячеек лучи используются как средство достижения точности при расчете коэффициентов распределения энергии точечного источника, распределения — по ячейкам сети.
Для определения значения подынтегральной функции проводится геометрический расчет пересечения ячеек многими лучами. Это существенная особенность всего рассмотрения. Она возникает при переходе от расчета коэффициентов видимости — для переноса в вакуумном приближении — к расчету коэффициентов
ослабления — для переноса с учетом взаимодействия излучения и вещества.
Величину импульса можно вынести за знаки интегралов правой части интегрального уравнения и сократить на эту величину обе части уравнения. Остается интегральное описание распределения импульса излучения, равного единице. Справа при этом под интегралами — функции распределения полной вероятности — единицы, функция распределения проникающего излучения по граничной поверхности и функция распределения выведенной веществом части излучения внутри объема расчетной области.
Для выбранного варианта описания переноса (нестационар- ность, приближение серой материи, квазистационарность границ и коэффициентов взаимодействия излучения вещества, изотропность рассеяния, локальное равновесие собственного излучения и вещества, закон Ламберта для излучения с границы области) функция распределения энергии точечного источника выражена в аналитической форме.
Последующие дискретные аналоги формул распределения отличаются заменой интегралов суммированием. Можно понять и проверить правомерность действий с бесконечно малыми величинами в предшествовшем изложении. Основные дискретные формулы описания распределения — алгоритмическая реализация аналитических формул, подготовка к написанию программы на алгоритмическом языке для ЭВМ — скоро не будет требоваться. Основная трудность в описании геометрии расчета, в построении сети ячеек. Мы приложили усилия для описания широкого класса двумерных сеток с разделением геометрической и топологической информации.
До изобретения ЭВМ использовались таблицы основных аналитических функций, механические и электрические настольные счетные машины. С изобретением ЭВМ появился круг квазиана- литических математических средств вычисления, использование которых выполнимо на ЭВМ с гарантированной заказанной точностью, [9]. К таким средствам вычисления в [9], как примеры, отнесены используемые нами итерации Ньютона для решения нелинейных алгебраических уравнений (метод касательной) и вычисление однократных определенных интегралов. В [9] речь шла о первых ЭВМ середины XX века. Мы используем итерации Ньютона на стадии замыкания вычислений для одного шага по времени — после описания распределения энергии сеточных точечных источников. Итерации Ньютона используются локально, для каждой ячейки отдельно, — при решении нелинейных
алгебраических уравнений граничных условий диффузионного типа [7] и уравнения локального баланса вещества [8]. Мы используем приближенное вычисление экспоненты и двукратные определенные интегралы — по поверхности, а точнее — по телесному углу.
Когда-то в XIX веке встал вопрос о разработке экспериментальных установок с разрешающей способностью для установления 6 точных значащих цифр при измерении результатов опыта. Началось с опытов А. Майкельсона по измерению скорости света. В XX веке был изобретен, например, электронный микроскоп.
При современном прогрессе в развитии вычислительной техники необходима и возможна разработка численных моделей описания энергопереноса излучением на уровне средств геометрической оптики, с использованием имеющихся аналитических решений, таких как закон Бугера-Ламберта при описании переноса вдоль луча, — с гарантированным определением 6 и более точных значащих цифр для результата численного эксперимента.
На первых тестовых расчетах с помощью первых программ, основанных на изложенном новом подходе, достигалась точность результата в 6 и более значащих цифр.
Далее более подробно будет описана функция углового распределения импульса излучения точечного источника для фиксированного луча и ослабления на пути источник-цель. Это алгебраические формулы с использованием функций косинуса и экспоненты.
Для описания переноса рассматривается пучок расходящихся лучей с центром в точке-источнике Р, точнее рассматриваются лишь лучи этого пучка, выводящие внутрь расчетной области. Говорят об актуальной части ^(Р) телесного угла, связанной с фиксированным точечным источником Р. Для внутреннего точечного источника ^(Р) = 4л. Для граничного точечного источника ^(Р) = 2л, если точка принадлежит гладкому участку граничной поверхности. В противном случае, если граничный точечный источник принадлежит линии излома граничной поверхности, величина актуальной части ^(Р) телесного угла будет иной. Мы договорились о том, что точечных источников, расположенных на линии излома, рассматривать не будем.
Рассмотрим геометрические лучи PQ, PR; PQ 6 ^(Р), Q Е S, PR Е ^(Р), R Е V. Говорят, что точка-источник Р и граничная точка-цель Q видят друг друга, если на отрезке PQ от
источника до границы нет других точек пересечения с границей, кроме точки Q; говорят также, что точка Q принадлежит области видимости S^P) на граничной поверхности S для точечного источника Р, Q Е S(P), S(P) С S. Говорят, что точка-источник Р и внутренняя точка-цель R видят друг друга, если на отрезке PR от источника до внутренней точки нет точек пересечения с границей; говорят также, что точка R принадлежит области видимости V(P) внутри объема V для точечного источника Р, ReV(PY V(P) С V.
Для фиксированного точечного источника Р выполняется расчет распределения импульса излученной им энергии:
— расчет углового распределения по направлениям расходящихся лучей описания переноса — геометрических лучей PQ с учетом приписанных к ним окрестностей телесного угла, или
JPQPQ1
этот расчет связан с вычислением областей видимости S^P\ V^P); вводится понятие — луч описания переноса — совместно геометрический луч PQ и его окрестность телесного угла, du:^ или Aw-^;
расчет распределения энергии при прохождении каждого фиксированного луча описания переноса: расчет распределения энергии, выведенной веществом, при пересечении внутренних ячеек; расчет распределения в граничных ячейкахпроникающей части импульса энергии;
лучи PQ, пересекающие границу S до прихода в точку Q, выводятся из рассмотрения; в этом состоит учет затенения.
В случае невыпуклости граничной поверхности для фиксированной точки-источника может произойти затенение отдельных участков граничной поверхности и объема расчетной области. Внутри актуальной части телесного угла (2тг или 4л) исключаются из расчета лучи, описания переноса, ведущие в затененные участки границы. После первого пересечения геометрического луча PQ с границей S и выхода его за границы расчетной области V описание переноса вдоль этого луча прекращается, он выбрасывается из рассмотрения.
Фиксируем теперь и точку-источник Р, и луч PQ. Будем называть его актуальным лучом, если он ведет в точку-цель Q области видимости S^P) на границе, Q Е S^P). Каждому актуальному лучу PQ ставим в соответствие некую окрестность
телесного угла: duu^ — при аналитическом описании, ^шрд — при дискретном описании.
Для фиксированной точки-источника Р будем называть актуальный луч PQ и окрестность телесного угла, Лш^ или Ди—, лучом описания переноса.
Надо добиться, чтобы сумма дискретных долей ^^р^ Для совокупности всех дискретных лучей PQ описания переноса для фиксированной точки-источника Р равнялась актуальной части Q(F) телесного угла,
для граничного источника Р —
£Д^ = 2тг,(26)
PQ
для внутреннего источника Р-
£Д^ = 4тг.(27)
PQ
Мы достигнем этой цели, с одной поправкой, если при разбиении телесного угла Q(P) будем использовать элементы разбиения граничной поверхности, о поправке ниже. Производится разбиение граничной поверхности S на элементы AS, S = идз. Геометрическим местом размещения точек-целей Q избираем множество счетных центров Qas элементов AS. Каждому лучу PQas ставим в соответствие в качестве его окрестности телесный угол, определяемый конической поверхностью с центром в точке Р, опирающейся на границу элемента разбиения AS со счетным центром Qas-
Поправка такова. В случае невыпуклости поверхности S и возникновения затенения нужно разбивать на элементы AS не граничную поверхность S, а область видимости 8(Р\ Это дорого, но возможно при решении двумерных задач с осевой симметрией с предварительным использованием алгоритма, описанного в [16], [21], аналитического расчета дуг видимости. Использование этого алгоритма приводит к уменьшению на порядок величины ошибок вблизи границ области видимости. Без использования этого алгоритма на границах области видимости величина коэффициентов углового распределения завышается на величину первого порядка малости по сравнению с половиной площади элемента AS разбиения. Ошибки контролируются и исправляются делением коэффициентов распределения на их
сумму. Сумма исправленных (нормированных) коэффициентов распределения равна единице.
Обсудим возможности применения для разбиения телесного угла двух формул для определения элемента телесного угла deep, являющегося окрестность луча PQ. Первая формула дает описание пучка расходящихся лучей PQ с центром в точке Р, используя местную сферическую систему координат,
dup = smddddcp. (28)
Для граничной точки-источника Р за ось сферической системы координат удобно выбрать вектор пр нормали к граничной поверхности S, внутренней по отношению к области V расчета. Формула (28) определяет лучи описания переноса локально для Р, исходя из информации лишь о точке-источнике. Лучи PQ никак не связаны с точками-целями на граничной поверхности S. Встанет вопрос о распределении части энергии, проникающей вдоль луча описания переноса, между граничными ячейками, соседними с точкой Q, Q Е S(P\ Ведь пересечение границы S с конусом телесного угла dwp, в общем случае, не совпадает с площадью граничных ячеек.
Подобным образом определяют лучи описания переноса в разностных схемах, основанных на аппроксимации интегро-дифференциального уравнения переноса.
Вторая формула дает описание пучка расходящихся лучей PQ с центром в точке Р, используя проекции граничных ячеек на поверхность сферы единичного радиуса,
duP = -^-dSQ, (iQ = cos (uq^QP). (29) Ppq
Напомним, что телесный угол измеряется отсеченной им площадью поверхности сферы единичного радиуса, подобно тому, как линейный угол меряется отсеченной им длиной дуги окружности единичного радиуса. Здесь вектор пу нормали к граничной поверхности S, внутренней по отношению к области V, рассматривается для точки Q границы, для элемента dSq площади границы и одновременно для точки Р источника. Дискретный вариант формулы (29) для телесного угла Awp, опирающегося на граничную ячейку Sq,
^шр = . (30)
Ppq
Можно рассматривать дополнительное разбиение ^SqiSq = = |J /\Sq граничной ячейки Sq и численно вычислять интегралы по телесному углу, опирающемуся на граничную ячейку. Этим мы будем пользоваться при вычислении с помощью (30) коэффициентов видимости, например интеграла
[ ^Sq.
Sq(P) PPQ
Цель достигнута.
Остановимся на формуле (28), которую можно воспринимать как элемент площади на сфере единичного радиуса — элемент площади поверхности вращения радиуса sin?? — вокруг оси с направлением пр внутренней нормали к касательной плоскости граничной поверхности S в точке Р.
Угловое распределение одностороннего потока dEP^ = J+(P,t)dSPdt
для граничного источника Р по направлениям PQ описывается стандартной формулой для потока векторной величины, только интегрирование ведется в пределах половины телесного угла,
j+(P3)= j I^P,PQ,t^ cos(2tt)
Формула (28) позволяет аналитически описать угловое распределение одностороннего потока — для поверхности, излучающей по закону Ламберта, т. е. с изотропной интенсивностью I, не зависящей от направления PQ,
I(P,PQ,t) = I(P,t) = Р
Для интегрирования по полусфере телесного угла в точке-источнике Р границы мы рассмотрим ту же сферическую систему координат (с полярной осью в направлении пр), какую мы рассматривали в связи с (28),
2тг тг/2
J+ = I Г Г cos 1? sin i? da) dtp = Itt, I = — J+,J J7Г
0 0
— с заменой переменных cost? = /z, sint? = —dpt.
В случае, когда поверхность излучает по закону Ламберта, интенсивность в направлении FQ, для которого угол наклона вектора PQ к вектору нормали пр равен #,
= Icos^ = - J+(F,F.
7Г
В качестве упражнения предлагаем читателю тем же способом найти нормирующий множитель при описании углового распределения излучения граничного точечного источника в случае усиления преимущества направлений PQ, близких к нормали пр. Пусть функцией распределения граничного излучения является не cos??, как для поверхностей, излучающих по закону Ламберта, a cos3 тР
В практике счета трехмерных задач переноса излучения внутри полостей встречались случаи описания распределения в виде таблицы. Указывались диапазоны изменения cos?? и соответствующие им значения весовой функции. В алгоритме предусмотрена нормировка коэффициентов делением на их сумму.
Перейдем к сравнению основных особенностей подхода к численному решению краевой задачи на основе фотометрической модели и на основе интегро-дифференциального уравнения и интегрального уравнения Пайерлса.
Правая часть интегро-дифференциального уравнения описывает точечный источник Р с помощью сложного функционала: для излучаемой энергии, интегральной по направлениям, надо указать часть, излучаемую в фиксированном направлении. При этом размерность правой части уравнения должна совпасть с размерностью левой части, описывающей изменение интенсивности в процессе переноса.
При численной аппроксимации интегро-дифференциального уравнения на сети ячеек возникают трудности, связанные со сложностью аналитического языка описания физического процесса переноса: трудности согласования счетных неподвижных направлений переноса с сетью движущихся счетных ячеек, трудности согласования описания источников излучения и переносимых долей энергии. Используется трудно контролируемый математический аппарат совместного решения всех аппроксимирующих уравнений с дополнительными итерациями этого решения на каждом шаге по времени. Фактически дискретное описание переноса на основе интегро-дифференциального уравнения не есть следствие аналитических средств, предоставляемых самим уравнением. Привлекается опыт решения другого уравнения —
уравнения лучистой теплопроводности в диффузионном приближении. Вводится в рассмотрение величины потоков энергии через грани ячеек. Привлекается так называемое «дополнительное соотношение», означающее интерполяцию по сдвиговым переменным и выход из класса гладкости решения уравнения переноса.
Интегральное уравнение Пайерлса и интегро-дифференциальное уравнение эквивалентны, они описывают один и тот же физический процесс переноса излучения в одних и тех же предположениях. При использовании нестационарного уравнения Пайерлса делается дополнительное предположение о квазистационарности геометрии расчетной области, об этом предположении — чуть ниже. Аналитически уравнение Пайерлса одновременно и проще, и богаче. Проще потому, что рассчитываемые функции имеют меньшую размерность, чем интенсивность: они интегральны по направлениям переноса, зависят от меньшего числа переменных. Богаче потому, что описывается зависимость плотности потока энергии, приходящей в точки-цели, от плотности потока энергии, излучаемой всеми точками-источниками с учетом геометрии расчетной области, с учетом запаздывания и затенения. Несмотря на все это при численном описании переноса на основе уравнения Пайерлса не удается избежать совместного решения всех аппроксимирующих уравнений с дополнительными итерациями этого решения на каждом шаге по времени. Само уравнение расщепляет вычисления по точкам- целям, но не расщепляет их по точкам-источникам.
Для граничных и внутренних точек-целей описывается перенос к этим точкам проникающей части излучения от всех точек- источников, граничных и внутренних. Не доходит та часть первоначального излучения источников, которая выведена веществом за счет поглощения и рассеяния. В свою очередь, точки-источники излучают лишь совокупную энергию фотонов, «родившихся» в окрестности этих точек-источников — за счет рассеяния и за счет собственного излучения вещества в предположении локального равновесия между поглощением и излучением (закон Кирхгофа).
В первой программной реализации нового метода [33-35] и в данной публикации используется предположение о квазистационарности границ счетных ячеек и величин коэффициентов взаимодействия излучения и вещества при счете для одного шага по времени [7]. Предположение о квазистационарности границ расчетной области всегда использовалось и ранее в нестационарных расчетах — при вычислении угловых коэффициентов [17, 24]
в расчетах переноса излучения для полостей сложной формы на основе уравнения Пайерлса [12, 15, 17]. Учитывалось, что положение границ расчетной области мало меняется за время прохождения света от источника до цели вследствие большого отличия скорости света и скорости вещества.
Предположение о квазистационарности границ обеспечивает квазистационарность величин коэффициентов рассеяния, зависящих от плотности вещества. Ведь из квазистационарности границ следует квазистационарность плотности вещества.
Но из квазистационарности границ следует и квазистационарность температуры вещества. Предположим от противного, что температура вещества меняется со скоростью переноса проникающего излучения. Тогда с такой же скоростью должно меняться давление вещества, ускорение, скорость и границы вещества.
Упомянутые предположения квазистационарности дают возможность рассматривать интегралы по поверхности и объему, определять область видимости для точечного источника, описывать распределение импульса фиксированного источника на всю глубину области видимости.
Опираясь на предположения квазистационарности, мы будем опускать аргумент времени для функций состояния вещества, коэффициентов взаимодействия излучения и вещества и для функции оптической толщины, считая, что время для этих функций совпадает со временем описания источников,
ПРЛ\ р^РЛ\ WHY
Xa(P,t\ X^P^Y Xs(P,t), Xa,s^P,t\ Xr,s^P,t\ TpQ(f).
К этим функциям не будут применяться сдвиги по времени, связанные с временем прохождения со скоростью света некоторого расстояния, например от источника до цели. В дискретном варианте описания на каждом шаге по времени используются значения этих функций, вычисленные на предыдущем шаге по времени.
Найдена новая интегральная форма уравнения корпускулярного переноса. Ее использование облегчает численное решение нестационарной краевой задачи видимости, позволяет вести счет с большей точностью и с большими возможностями детального контроля точности. Описание переноса носит (явный) характер последовательного распределения величины излучаемой энергии во времени и в пространстве: по времени излучения, по источникам, по лучам переноса, по пересекаемым элементам объема для
выведенной веществом части энергии, по элементам граничной поверхности для проникающей части переносимой энергии, по времени запаздывания. На каждом этапе распределения имеется возможность контроля локальной точности распределения и возможность соблюдения абсолютной балансности счета с помощью нормирования дискретной функции распределения энергии. Описание переноса сводится к умножению на числа, в сумме равные единице, к сложению порций энергии, одновременно перенесенных от разных источников. Все эти арифметические операции выполняются со всеми вычисляемыми ЭВМ точными значащими цифрами
При замыкании вычислений в каждой внутренней ячейке локально итерациями Ньютона решается алгебраическое уравнение четвертой степени.
При счете коэффициентов со всеми вычисляемыми ЭВМ точными значащими цифрами решается геометрическая задача пересечения луча и поверхностей вращения второго порядка. Так рассчитывается подынтегральная функция при численном вычислении определенного интеграла — по телесному углу.
Сначала на основе метода коэффициентов видимости удалось построить явную автоматически балансную схему [18] численного моделирования нестационарного лучистого энергопереноса внутри полостей сложной формы с неподвижными границами в вакуумном приближении, без учета взаимодействия излучения и вещества. Это дало возможность с 1972 года перейти к счету трехмерных задач этого типа. С 6 точными значащими цифрами в разнообразных вариантах постановки граничных условий решается, в частности, трехмерная стационарная тестовая задача В. Г. Рогачева о переносе излучения через сферическую полость, имеющая аналитическое решение. По программе [18] было решено до тысячи трехмерных задач.
Затем, на основе алгебраических законов фотометрической теории и теории локально равновесного переноса в приближении серой материи удалось построить алгоритмы численного решения двумерных нестационарных задач с учетом изотропного рассеяния, поглощения и локально равновесного собственного излучения вещества, [33].
При тестовой отладке модулей новой программы использовались сравнительные расчеты серии задач, допускающих ква- зианалитический счет с помощью пакета программ MathLab. Результаты расчетов совпадали в пределах задания точности для MathLab, до 6-8 точных значащих цифр, [37, 38].
Со всеми опубликованными 6 точными значащими цифрами на существенно двумерной счетной сети ячеек была решена известная тестовая одномерная нестационарная задача о выходе импульса излучения из центра шара с учетом изотропного рассеяния [10].
И, наконец, предельным переходом была получена новая форма аппроксимируемого интегрального уравнения переноса излучения, [34, 35]. Вместе с известным интегральным уравнением Пайерлса [2, 3] новое уравнение замыкает интегральное описание переноса через область на основе рассмотрения пучков лучей.
Алгоритм метода [33] глубоко расщепляется и удобен для использования параллельных вычислений.
Новый метод по характеру вычислений имеет много общего с методом Монте-Карло и сможет, по-видимому, совмещаться с ним в рамках единой программы.
Глава 1
ИНТЕГРО-ДИФФЕРЕНЦИАЛЬНОЕ УРАВНЕНИЕПЕРЕНОСА ИЗЛУЧЕНИЯ
Уравнение в частных производных
Ограничим рассмотрение корпускулярного переноса излучения приближением серой материи, учетом изотропного рассеяния, поглощения и собственного излучения в условиях локального термодинамического равновесия и локального равновесия между процессами поглощения и собственного излучения вещества.
Интегро-дифференциальное уравнение переноса излучения получено как односкоростной линеаризованный вариант кинетических уравнений Больцмана в предположении некогерентности источников. В обозначениях, близких к [7, 11], интегро-дифференциальное уравнение имеет вид:
D др, оз) = -pnr.s(P,£) + /(P,£)),(31)
D = -^ + (fi,grad) + xa,$.(32)
с ot
Нижние индексы в (31), (32) напоминают о физической природе коэффициентов взаимодействия xQvS и функции источника в правой части уравнения,
(33)(33)хт^Р, i) = х^Р, i) + XS(P, i),nr,s(P, t) = uT(P, i) + us(P, t\ur(P,t) = хасаТ4(РД\us(P,t) = xscU(P,t),cU(P,t) = | I(P^t)duj.
(33)
(33)
(4тг)
Здесь I(P, fi,t) — интегральная по частоте интенсивность — скалярная функция, описывающая энергоперенос излучением
в момент времени t в направлении FQ луча, выходящего из точки-источника Р. Даже после интегрирования по частоте интенсивность является функцией шести независимых переменных: трех координат точки-источника Р, двух угловых параметров, определяющих направление луча FQ, и времени t.
Функция интенсивности позволяет описать угловое распределение плотности потока энергии точечного источника: внутри элемента duo телесного угла вокруг фиксированного луча PPL эта плотность равна ЦР^Д^Дси.
Размерность интенсивности: [/(Р, fM)] = Вт/(м2сср), это плотность направленного потока энергии — величина энергии, переносимой в направлении вектора РР1 в единичном телесном угле (ср — стерадиан) вокруг этого вектора, через единицу площади, в единицу времени.
Размерность дифференциального оператора [D] = 1/м.
Левая часть уравнения (31), (32) содержит дифференциальный оператор в частных производных. Описывается изменение интенсивности излучения при переносе вдоль фиксированного луча.
Член уравнения, содержащий коэффициент полного взаимодействия xa,s, xaiS = ха + xs, описывает ослабление излучения веществом с учетом поглощения и рассеяния. Нижний индекс а указывает на связь величины с поглощением, нижний индекса — с рассеянием. В условиях локального термодинамического равновесия коэффициенты взаимодействия ха, xs, ха^8 излучения и вещества — заданные функции состояния вещества. Размерность \иа,И8,ИаД = 1/м.
Размерность членов уравнения (31), (32) и его правой части Вт/(м3сср). Это размерность направленной мощности источника, размерность энергии, излучаемой в единицу телесного угла в окрестности направления PPL единицей объема вещества в единицу времени.
Правая часть уравнения (31), (32), интегральная по направлениям с учетом (33), описывает изотропную функцию точечного излучения и вклад этого точечного источника в направление PPL за счет собственного излучения вещества (нижний индекс г), за счет изотропного рассеяния (нижний индекс s) и с помощью заданной функции излучения f(P,t). Последний вклад, f(P,t\ будем в дальнейшем изложении опускать для простоты написания, считая, что функция источника / добавляется к описанию источника излучения вещества, если она не ноль. Величина плотности энергии излучения Ц\РД) пропорциональна интегралу интенсивности по всем направлениям переноса.
Правая часть уравнения (31), (32) описывает угловое распределение энергии для точечного источника излучения Р.
Размерность произведения cU(P,t) — размерность плотности потока энергии в окружающее пространство — [ct7(P, £)] = Вт/(м2с), с — скорость света. Размерность выражения —- U(P,P совпадает с размерностью интенсивности.
Перейдем от уравнения (31), (32) к дифференциалу, тогда xa)Scdt — доля излучения, выведенного веществом за время dt. Соответственно, xacdt — доля излучения, поглощенного веществом, xscdt — доля излучения, рассеянного веществом, — за то же время.
В случае изотропного рассеяния — xsU(P,t) — энергия излучения, рассеянная единицей объема вещества в единицу времени внутрь единичного телесного угла.
Введем обозначение
us(P,t) = cxsU(P,t\(34)
— мощность точечного источника Р за счет рассеяния излучения веществом — количество энергии, рассеиваемой в окружающее пространство единицей объема вещества за единицу времени, размерность [ns] = Вт/(м3с).
Использование формулы (34), в частности, приводит к численным схемам, построенным на аппроксимации интегро-дифференциального уравнения, на совместном решении системы всех аппроксимирующих приближенных уравнений и к итерациям для определения функции полной плотности U(P,tY без чего нельзя описать точечные источники, и§1 иг, для счета следующего шага по времени.
Независимыми термодинамическими параметрами вещества (для газов) будем полагать плотность р и температуру Т вещества внутренних ячеек. Будем считать значения этих величин функциями пространственных координат и времени, P(p,tCT(p,ty
Будем считать заданными функциями температуры и плотности вещества:
коэффициент поглощения и коэффициент рассеяния
ха = Ш,р\ xs = fM(35)
плотность внутренней энергии вещества (нижний индекс е)
Ue = CvpT;(36)
размерность плотности внутренней энергии вещества Ре\ = Вт/М3,
задан удельный коэффициент Су теплоемкости вещества;
мощность точечного источника (Р, t) за счет собственного излучения вещества;
количество энергии, испускаемой самопроизвольно в окружающее пространство единицей объема вещества за единицу времени,
иТ = схааТ\(37)
здесь аТ^ — плотность локально равновесного излучения абсолютно черного тела (закон Стефана-Больцмана), постоянная Стефанаа задана; коэффициент поглощения ха играет роль коэффициента «серости» вещества (закон Кирхгофа). Размерность плотности собственного излучения вещества \ит^ = Вт/(м3с).
Правая часть уравнения описывает функцию источника — излучательную способность точечного нестационарного источника — и связь этого источника с направлением FQ. Интегро- дифференциальное уравнение имеет вид:
D ЦР, Ci,t) = ф ит ЦР,(38)
Характеристическое уравнение
Дифференциальная часть оператора D, (32), (38) — развернутая форма полной «субстанциональной» производной — по длине g пути излучения вдоль луча FQ от точки Р, dg = cdt, переменная времени пропорциональна пути. После перехода
к характеристической переменной g интегро-дифференциальное уравнение (31), (32) имеет вид обыкновенного дифференциального уравнения первой степени:
[^ + Xa's^\ Z^9’Р^' t + g/c) = ^ urAp^ Г (39)
Уравнение (39) — эквивалентная форма интегро-дифференциального уравнения (31), уточняющая его. В характеристической форме (39) дифференциальный оператор D интегро-дифференциального уравнения (31), (32) совпадает с дифференциальным оператором одномерного уравнения Ламберта, см. ниже (43).
Уравнение (39) наиболее удобно для установления класса гладкости решения: интенсивность дифференцируема в направ-
лении луча, выводящего из точечного источника. Нет непрерывности по сдвиговым переменным и по углам.
Хорошо известно, что такое выглянуть из-за угла. Можно увидеть что-то непредсказуемое. Так же и видящий бег по своему телу пятнышка от наводящего прицел лазера не знает, где ствол, если не доведется взглянуть вдоль лазерного луча.
Алгебраическое уравнение
Мы переходим к третьей, эквивалентной и самой ранней по времени открытия, форме интегро-дифференциального уравнения переноса, — к алгебраической форме закона Бугера-Ламберта.
Обычно решение обыкновенного линейного дифференциального уравнения (39), и см. ниже (43), (44), записывают в виде:
MQ, PQ, t) = exp (-tpqW, PQ, t - pPQ/Q. (40)
Формула (40) читается так: интенсивность I(Q,PQ,t) в точке Q в момент времени t по направлению PQ включает в себя экспоненциально ослабленную интенсивность излучения, рожденного в предшествующий момент времени, с учетом запаздывания, в точке-источнике Р с этим же направлением. Нижним индексом Р мы сообщаем место рождения части излучения. В интенсивность I(P, PQ,t — ppq/c), кроме интенсивности, обеспеченной точечным источником ur^P,t — ppq(c\ может входить на правах слагаемого ослабленная интенсивность, описываемая аналогично (40), — рожденная в том же направлении предстоящими Р точками R — в еще более ранние моменты времени t — Ppq/c — Prp/c, Р Е RQ, RP + PQ = RQ. В интенсивность I(Q, PQ, t) могла бы входить на правах слагаемого и ослабленная интенсивность из последующих Р точек Rf,Rf Е PQ.
В формуле (40) нижний индекс Р указывает на вклад в интенсивность только точечного источника Р.
Экспонента (безразмерная величина) решение однородного уравнения
{4 + ха>Лр\ех.р [-т(Л] = 0,
д — длина пути излучения вдоль луча РП от точки Р, dg = cdt.
Степень tqp экспоненты (безразмерная величина) оптическая толщина пути QP, измеренная в пробегах. Пробегом называ-
ется путь излучения, при прохождении которого интенсивность уменьшается в е раз.
gTQg^ =QgTQg^ =Q7-dg = [ xaiS(g)dg,Q7-dg = [ xaiS(g)dg,QПробег la,s и коэффициент взаимодействия xa)S излучения и вещества — взаимно обратные величины, la,s = l/xa,s,
g
TQg^ =
Q
g
TQg^ =
Q
7-dg = [ xaiS(g)dg,Q
7-dg = [ xaiS(g)dg,Q
у [ ^a,s(9) dg = иа1Цд\
9 Q
g — координата длины пути,
TQP
Размерности: [y,Z»,s] = м, \xaiS\ = 1/м.
Вообще говоря, интенсивность Ip(Q, PQ,t), обязанная своим рождением точечному источнику Р, может быть лишь слагаемой интенсивности I{Q, PQ,t).
Нередко в изложении остается опущенным одно из двух предположений исходной экспериментально-теоретической корпускулярной фотометрической модели переноса излучения о взаимном влиянии точечных источников друг на друга:
интенсивности, рожденные разными источниками с разными направлениями, не влияют друг на друга (вспомним пример опыта с камерой-обскурой),
интенсивности, рожденные разными источниками с одинаковыми направлениями, складываются. Днем выключенная фара видна за счет солнечного излучения, включенная днем фара видна с суммарной интенсивностью.
Решение уравнения (39) мы будем использовать в одной из двух форм, с учетом нестационарности задачи:
Ip(Q,PQ,t) = exp(-TPQyriS(P,PQ,t - ppq/c), (41)
Ip(Q, PQ, t + ppq/c) = exp (-TPQ)Jr,s(P, PQ, t\ (42)
В отличие от (40), формулы (41) и (42) подчеркнуто разделяют интенсивность по «генетическому» признаку, описывая перенос квантов от места их рождения.
Формулы (41) и (42) отличаются заменой переменной времени.
В (41) фиксированный момент времени t приписывается точке-цели Q — при описании переноса вдоль луча PQ излучения,
рожденного точечным источником Р в предшествовавший момент времени t — Ppq/c. Описывается влияние на интенсивность в фиксированной точке-цели Q излучения, рожденного в направлении PQ точечным источником Р.
В (42) фиксированный момент времени t приписывается точке Р — точечному источнику энергии излучения. В результате переноса вдоль луча PQ в последующие моменты времени t + Ppq/c энергия, рожденная этим источником в момент i, дойдет до точки-цели Q экспоненциально ослабленной. Описывается вклад излучения, рожденного фиксированной точкой- источником Р, на интенсивность излучения в точках-целях Q для фиксированного луча PQ.
Экспоненциальный закон ослабления интенсивности за счет взаимодействия излучения и вещества среды переноса был установлен опытным путем французским ученым П. Бугером в 1729 г. и был получен немецким ученым И. Г. Ламбертом в 1760 г. на основе теоретической модели, обобщающей экспериментальные результаты. Закон называется законом Бугера- Ламберта [1, 22]. Формула (40) была получена Ламбертом как решение одномерного уравнения для луча, эквивалентного уравнению переноса (39) в характеристической форме,
у = -xdg, (43) в нашем случае
X = ха,8. (44)
В форме (41) закон Бугера-Ламберта имеет привычный вид решения уравнения (39). В форме (42) тот же закон является третьей и наиболее ранней алгебраической формой уравнения переноса, эквивалентной интегро-дифференциальному уравнению переноса. Описывается распределение лучистой энергии от источника вдоль луча переноса с учетом запаздывания.
Формула (32) придает интегро-дифференциальному уравнению (31) вид уравнения в частных производных, содержащего производные первого порядка, формула (33) придает правой части уравнения интегральную форму. В (38) интегральной является величина плотности и8 рассеянного излучения (интеграл по телесному углу) — совокупная плотность фотонов, поменявших свое направление в единице объема вокруг точки Р в единичном диапазоне времени относительно фиксированного момента времени t.
Формула (39) уточняет формулу (32): любой луч пучка расходящихся лучей PQ является характеристикой уравнения (31). Уравнение (31) вдоль любой характеристики-луча становится обыкновенным дифференциальным уравнением первого порядка. Другого содержания, отличного от указанного, уравнение в частных производных (31), (32) не имеет. Отметим, что статья энциклопедии [27] приводит именно характеристическую форму (39) уравнения в качестве интегро-дифференциального уравнения переноса излучения.
Формулы (41), (42) описывают проникающую часть излучения точечного источника, достигающую фиксированной точки на любом луче.
Формула (42) дает описание переноса излучения от точечного источника вдоль луча, в еще более уточняющей третьей форме, алгебраической, позволяющей проводить непосредственные вычисления. Уравнение переноса в форме (42) можно характеризовать как интегро-алгебраическое уравнение. Левая часть уравнения описывает распределение вдоль луча энергии, излучаемой интегральным по направлениям точечным источником Р в фиксированном направлении луча PQ. Правая часть, содержащая интеграл по направлениям, описывает долю энергии, излученной в данном направлении PQ и ее экспоненциальное ослабление на пути к цели Q. В нашем случае рассмотрения изотропного излучения все направления равновероятны, нужно привести правую часть уравнения к единице телесного угла.
С помощью формулы (42) можно описать распределение вдоль луча PQ энергии излучения фиксированного точечного источника Р — проникающую на границу часть излученной энергии и распределение вдоль луча выведенной веществом энергии. Напомним, что интегральное уравнение Пайерлса описания распределения выведенной веществом энергии не дает.
До рассмотрения общего случая пересечения лучом ячеек произвольной внутренней сети, предположим, что мы имеем дело с плоско-параллельным одномерным переносом фотонов, описываемым интенсивностью в двух направлениях PQ и QP, ортогональных слоистой системе ячеек вещества (приближение вперед-назад, [7, гл. 2]). Точки Р и Q принадлежат граничным плоскостям.
Мы не будем заниматься всеми вопросами полной постановки задачи в приближении вперед-назад, например, вопросами углового усреднения интенсивности. В нашу задачу входит
познакомиться с переносом вдоль одного луча PQ на основе формулы (42).
Будем считать, что граничный точечный источник Р в момент времени t за время т в единственном направлении PQ, перпендикулярном плоской границе, с единичной площади границы, Sp = 1, излучает импульс Е^РД\ пропорциональный диапазону времени т и плотности J+(F,t) одностороннего потока энергии в направлении PQ, Е^РД) = J+(P,t)Spr.
Опишем с помощью (42) распределение импульса Е(РД) на луче PQ — внутри цилиндра с единичной площадью сечения — с выделением по пути сначала трех отрезков — слоистых ячеек описания вещества — PRx, R\R^, R^Q. Обозначим tprx, trxr2,tr2q оптическую толщину слоев.
Точка входа в первый слой — точка Р излучателя. Энергия на входе в первый слой — энергия Ер = Е(РД), введенная излучателем.
Точка выхода из первого слоя — точка Rx пересечения луча с последней границей этого слоя. Энергия на выходе из первого слоя, (42), с учетом оптической толщины trrx этого слоя — энергия Erv ERi = EPexp(-TpiP.
Энергия EaiSPRi, выведенная веществом первого слоя PRx, — разность введенной и вышедшей энергии для этого слоя,
Ea,sPRx k^E рх)^ = Ер - Erx , рх = 0,5 \Р Rx |. (45)
Точка входа в очередной слой — точка выхода из предыдущего слоя — точка Rx. Энергия на входе в очередной слой — энергия на выходе из предыдущего слоя — Erx.
Точка выхода из очередного слоя — точка R2 пересечения луча с последней границей этого слоя. Энергия на выходе из очередного слоя — с учетом оптической толщины trxr2 этого слоя — энергия Er2, Er2 = Erx ехр(-тд,д2).
Энергия EagRxR2, выведенная веществом очередного слоя RxR2^ — разность введенной и вышедшей энергии для этого слоя,
Ea,sRxR2(t + р2/с) = Erx - Er2, р2 = PRx + 0,5 |FiR2 I. (46)
Точка входа в последний приграничный слой — точка выхода из предыдущего слоя — точка R2. Энергия на входе в последний слой — энергия на выходе из предыдущего слоя — Er2.
Точка выхода из последнего слоя — граничная точка Q окончания луча, последняя граница слоя. Энергия на выходе
из текущего слоя — с учетом оптической толщины tr2q этого слоя — энергия Eq, Eq = Er2 exp Ptr2qY
Энергия Ea)Sp2Q, выведенная веществом последнего слоя RsQ, — разность введенной и вышедшей энергии для этого слоя,
Ea,sR2Q^ + Рз/с) = Ер, - Eq, рз = Р R^ + 0,5 |Д2 Q\ • (47)
В результате переноса (45)-(47) энергия предстает в виде суммы четырех слагаемых,
E(J\t) = ^EajS{R, tp) Р Eq^P prq)c\ (48) R
Eas^R, tp) = Ea,sPRx + pi/c) + Easpxp2/t + P2Ic) +
R
+ Ea)SR2Q(t + рз/c),(49)
где слагаемое Eq — проникающая на границу часть введенной фиксированным источником Р энергии, остальные слагаемые описывают распределение внутри расчетной области V выведенной веществом части энергии.
Подставляя вместо величин энергии, выведенной слоями вещества, определяющие их разности, (45)-(47), получаем вместо (49) и (48):
У2 Е^№ tp) = Ер - Ерх + Ерх - Ер2 + Ер2 - Eq = Ер - Eq,
R(50)
E(P,t) = ЕР — Eq + Eq.(51)
Мы видим, (50), что в сумме у^ Ea)S(R, tp), определяющей R
распределение выведенной веществом энергии источника через разности входящих и выходящих потоков на границах каждого слоя, кроме первого и последнего члена правой части, все остальные входят в сумму и со знаком (+), и со знаком (-), и сокращаются. Эти члены описывают величины энергии, вышедшей из каждого слоя и, стало быть, введенной в следующий слой (для последнего слоя — на границу). Равенство (51) выполняется на ЭВМ с точностью вычисления экспоненты. Все сказанное о свойствах алгоритма описания переноса вдоль луча справедливо и при большем числе п пересекаемых слоев.
При п сю вместо (50)-(51) мы получим одномерное интегральное описание распределения вдоль луча PQ — внутри
цилиндра с единичной площадью сечения — импульса dE(P,t) энергии излучения граничного точечного источника F,
dE(P, t) = J+ (Р, t)SPdt, SP=i,
Q
dE(P,t) = j UaySp{g, t + g/c)dgdt + JpQ,t P ppq/c) dSq dt, p(52)
ua,spdg = dE^exp [~r(g - 0,5dy)] - exp [-r(y + 0,5^)]},
ea,spdg = -d{exp [-т(Л]}.(53)
Здесь J+(P,t), J~(Q,t) — плотности односторонних граничных потоков в направлениях PQ, PQ, соответственно; их размерность \J+(P,t), J~(Q, i)] = Вт/(м2с); I(g,PQ,t + д/с) = = dE exp [—т(д)] — решение характеристического уравнения, см. (43), (44),
\/fg + Xa’s^ ^g^^QHP д/с) = dE(P,t), экспонента ехр [—т(д)] — решение однородного уравнения; ua,sp(gP + д/с) — вклад точечного источника Р в линейную плотность поглощенной и рассеянной веществом энергии излучения; иа,8р(дЛ + g/c)dg — энергия излучения, выведенная веществом слоя dg, — доля импульса энергии dE(P,t) точечного источника Р; размерность [na,sp] = Вт/(м3с) — с учетом единичной площади сечения цилиндра; Jp(Q,t + ppq/c) dSQ dt — вклад точечного источника Р в линейный односторонний поток проникающего на границу излучения, доля импульса dE(P,t). Мы перешли к характеристической координате д вдоль луча PQ, отсчитывая ее от точки Р источника, луч PQ — характеристика, (39), интегро-дифференциального уравнения (31), (32). Выражение справа в (53) — та же величина, определенная на основе закона Бугера-Ламберта. Выражение внутри фигурных скобок и dg не зависят от энергии источника, поэтому элемент плотности энергии иа,8р(д, ^ + д/с) пропорционален излученному импульсу dE(P,t), точнее плотности одностороннего потока Л^^РД).
Функция еа18р(дД + д/с) — функция распределения энергии точечного источника Р, выведенной слоями вещества — линейная плотность поглощения и рассеяния излучения, соответствующая распределению единичного импульса.
Ниже, (55), мы будем рассматривать распределение вероятности быть поглощенной или рассеянной веществом и вероятности дойти до границы — для энергии излучения точечного источника Р. Единичный импульс — полная вероятность быть излученной для энергии точечного источника Р. Для импульса, равного единице, из (53) имеем:
eaAg, t + д/д = -"I" ехр[-т(д)].(54)
Для функций, связанных с описанием переноса излучения вдоль луча PQ, все аргументы пропорциональны пути д от источника Р, t = д/с.
QdE(P,t) = dE(P,t)рQdE(P,t) = dE(P,t)рИз (52) и (54) имеем:
Q
dE(P,t) = dE(P,t)
р
Q
dE(P,t) = dE(P,t)
р
dexp\ т(д^ ^ + dE(p, ^ ехр (_TpQ);
Q1 = I рQ1 = I рсокращая все члены на dE, получаем описание функции распределения граничного излучения вдоль луча:
Q
1 = I р
Q
1 = I р
(55)(55)dexp[-r(g)] , , z ч vv WJ dg + exp Ptpq^ .
(55)
(55)
^9
Единица — полная вероятность, вероятность излучить импульс энергии для точечного источника. Эта вероятность распределяется на две вероятности: вероятность поглотиться и рассеяться слоями вещества и вероятность достичь границы для проникающей части излучения.
Все части распределения полной вероятности в (55) выражены с помощью одного и того же аналитического средства — алгебраической формулы закона Бугера-Ламберта. При численных расчетах это дает возможность избежать дисбаланса энергии.
Пользуясь формулой интегрирования Ньютона-Лейбница, получаем подтверждение правильности (52), (55):
dE^P^t) =
= dE(P,t){exp(-TPp') - ехр (-тру)} + dE(P, Ц exp (~tPq), dE(P, Ц = dE(P, Ц — dE(P, Ц exp (—тру) + dE(P, Ц exp (—тру), 1 = exp (-tpp) - exp (-тру) 4- exp (-тру) = 1.
Во все слагаемые формул (48), (52) входит множитель At или, соответственно, т — шаг по времени. Если на него сократить, получатся формулы, описывающие распределение лучеиспускательной способности источника:
J+{P,t)Sp = ^Ua^R, tp^Vp + JP(Q, t + ppq/c)Sq, (56) R
Q
J+(P,t}dSp = j ua,s(g,t + g/c) dg + Jp(Q,t + ppq/c) dpQ. (57) p
Для сравнения с формулами (50)-(57) одномерного интегрального уравнения переноса излучения на основе формулы (42), приведем одномерное интегральное уравнение Пайерлса [2, 7] на основе формулы (41),
J (QH^Sq = J+(P,t - ppq/c)dSp exp(-TpQ) +
(58)(58)Q
(58)
(58)
2 «r,sO> t- д/с) ехр[-т(ж)] dg.
R
Здесь точка-цель Q в момент времени t собирает ослабленную согласно закону Бугера-Ламберта (41) энергию, высвечиваемую всеми источниками на луче PQ, граничным источником J+(P,t — ppq/c)dSp и внутренними источниками -и^дЛ — g/c)dg; - для приближения вперед-назад описывает
1
угловое распределение вместо —- для внутренних источников изотропного излучения.
В описаниях распределения на основе (42) один изолированный источник рассылает импульс своей энергии всем ячейкам вещества вдоль луча (доли поглощенной и рассеянной энергии) и граничной ячейке (доля проникающей части излучения):
Е(Р, ^=^ Е(Р, t^pp + Е(Р, t^pq,R
PR = exp (-TPR+) - exp (-tPr_), pQ = exp (-tPq\ (59) ^PR + PQ = УХехР(-Тиг+) -ехр(-трд_)} + exp (—тру) = 1, RR
(60)
Можно на основе (59) вычислить доли рр, pq распределения, проконтролировать точность вычисления суммы долей; в данном
случае точность будет на уровне точности счета ЭВМ, в результате чего счет будет балансным (60) — на том же уровне.
Для лучей, ведущих из остальных источников излучения, формулы описания переноса аналогичны.
В формуле (58) для уравнения Пайерлса нет аналитической основы для описания распределения выведенной веществом части энергии. Именно на этой основе замыкается расчет переноса излучаемой энергии — на шаге по времени.
Глава 2
ИНТЕГРАЛЬНОЕ УРАВНЕНИЕ ПАЙЕРЛСА
Это уравнение эквивалентно интегро-дифференциальному уравнению при решении краевой стационарной задачи внутри выпуклой области V, ограниченной замкнутой поверхностью S, [2, 3]. Вместе с тем, оно уточняет интегро-дифференциальное уравнение. Уравнения (31), (39), (42) описывали перенос вдоль фиксированного луча. Интегральное уравнение Пайерлса описывает перенос для пучка лучей, сходящихся в фиксированную точку-цель, с учетом рассматриваемой краевой задачи. Аналитически описывается функция распределения проникающей части энергии, переносимой вдоль лучей от различных источников, граничных и внутренних. При этом лучи пучка оказываются неравноценными при неизотропном угловом распределении. Примером служит излучение, описанное законом Ламберта.
Мы приведем интегральное уравнение Пайерлса, описывающее нестационарный перенос для краевой задачи внутри невыпуклой области с учетом видимости. Для точек границы описание переноса замыкается граничными условиями. Во внутренних точках совместно с уравнением переноса решается уравнение локального баланса энергии для вещества [8].
Мы не будем излагать какую-либо схему решения краевой задачи на основе уравнения Пайерлса. Пример такой схемы дается в [17]. Принципиально то, что использование уравнения Пайерлса не освобождает от решения полной системы приближенных уравнений и от выполнения глобальных итераций на каждом шаге для улучшения балансности счета.
Рассмотрим пучок лучей, сходящихся в точку-цель Q и заполняющих видимую из точки Q часть расчетной области K(Q), S(Q), не затененную невыпуклыми участками граничной поверхности S.
Из-за запаздывания точка Q связана с точками Р и Д, относящимися к разным моментам времени. Если рассматривать движущееся вещество, возникают проблемы с описанием геометрии
расчетной области, с учетом затенения и запаздывания. Поэтому при записи нестационарного интегрального уравнения переноса для области делают дополнительное предположение о квазистационарности геометрии расчетной области — предположение о стационарности геометрии на время переноса через область, т. е. на максимальное время описываемого запаздывания. Это дополнительное предположение опирается на большую разницу в величинах скорости света и скорости вещества.
Используя уравнение переноса вдоль луча (39), (41), напишем известное [2, 3] интегральное уравнение Пайерлса для области, связывая лучами точки-источники Р на граничной поверхности и точки-источники R внутри области вещества с фиксированной точкой-целью Q.
Уравнение (61) показывает его вывод — с учетом сопровождающих пояснений физического характера.
Для граничной точки-цели Q описывается перенос энергии от всех граничных источников излучения Р вдоль лучей PQ и от всех внутренних источников излучения R вдоль лучей RQ:
J"(Q,f) = | J+(P,t - pPQ/c)KpQexp(-TPQ)dSp +
S(Q)
+ | ur,s(R,t- pRQ/c)K,RQexp(-TRQ')dVR; (61) V(Q)
Kpq = ^^, pp = cos(nP,PQ), pQ = cos (nQ,QP), ^Ppq
Krq = , ^ ’ ^Q = cos (nQ’
Здесь J\ J- — плотности односторонних потоков лучистой энергии. Поток J+ описывает излучение граничных источников внутрь расчетной области. Поток J описывает излучение, приходящее на границу со стороны расчетной области. Верхний индекс (+) соответствует направлению вектора внутренней нормали к границе по отношению к расчетной области. Индекс (—) — обратному направлению. Размерность плотности односторонних потоков [J+, У-] = Вт/(м2с).
Предполагается, что односторонние потоки У+, J связаны между собой граничным условием, позволяющим определить энергию, высвечиваемую внутрь расчетной области, на основе
вычисленной энергии с противоположным направлением переноса [7].
Л' = !Р У
Обсудим физическое содержание аналитических выражений, составляющих подынтегральные функции.
К описанию источника относятся выражения
J+(P,t~ ppQ/c)dSp, иТДКД- pRQ/c)dVp.
Выражения Kpq ехр(—tpq), КгКдехр(—трд) описывают три фактора распределения энергии точечного источника излучения:
распределение по направлениям, PQ или RQ,
распределение внутри телесного угла по долям (окрест; ностям) телесного угла, отнесенным к направлениям, PQ или RQ,
ослабление на пути к точке-цели Q на границе.
Выражения Kpq, Krq служат для описания углового распределения излучения по направлениям PQ и RQ, соответственно; описываются первые два фактора распределения, распределение по направлениям и распределение внутри телесного угла.
Множитель pip в выражении Krq для излучения источника Р граничной поверхности определяет вероятность получить направление PQ: эта вероятность пропорциональна рр. Предполагается, что граничная поверхность излучает по закону Ламберта — пропорционально косинусу угла наклона направления к вектору внутренней нормали. Для изотропных внутренних источников R все направления равновероятны, отсутствует множитель, описывающий распределение излучения по направлениям RQ.
Совпадающая часть, (iq)Ppq, iAq! Prch выражений Kqp, Krp описывает распределение излучения по телесному углу. Определяется телесный угол, под которым из точки-источника, Р или R, видна единица площади граничной поверхности вблизи точки-цели Q. Единицы площади перекрываются, поэтому в расчетах углового распределения этой формулой не пользуются.
На основе этих формул нельзя сохранить баланс расхода энергии для сеточных источников при описании переноса в сторону разных ячеек-целей.
Если умножить обе части (61) на dSq, то выражения правой части уравнения dSg^tg/р2Рд и ASq^q/р^д определяют телесные углы, под которыми из точки-источника, Р или R, виден элемент
площади ASq граничной поверхности — окрестности точки-цели Q. Умножим еще обе части (61) на dt.
После этого уравнение Пайерлса описывает обмен импульсами между локальными источниками и локальной целью;
J+(P,t- ppQ/c)dSpdt, ur,s(R,t- Prq/c) dVRdt
— импульсы энергии, излучаемые точечными источниками Р и R; J~ (Q,t) dSQ dt — импульс энергии излучения, проникающего на границу в окрестности точки-цели Q. Выражения
ASqpqASqPq
2 И "л2
^Ppq^Prq
определяют угловые коэффициенты, называемые локальными, [11]. Множители 1/л, 1 /(4тг) — нормирующие множители; интегралы по всем актуальным направлениям от функций локальных угловых коэффициентов равны единице.
Третий фактор распределения энергии излучения точечного граничного источника — ослабление: до цели дойдет доля, ехр(—tqp) и ехр(—трр), первоначально излученной энергии — на основе закона Бугера-Ламберта. Доля, которая не дойдет до цели и будет выведена веществом,
l-exp(-TQp) и 1-ехр(-тДр).
Уравнение Пайерлса не описывает распределение вдоль луча выведенной веществом доли энергии излучения.
Уравнение (61) представляет собой интегральное уравнение Пайерлса переноса проникающей части энергии в граничные точки-цели Q.
Интегральное уравнение Пайерлса переноса энергии во внутренние точки-цели Q аналогично:
cu(Q,t) = | J+(P,t - ppQ/c]KpQexp-TFQ) dSp +
+ j ur,s(R,t- pRQ/c)KRQexp(-TRQ)dVR. (62)
KpQ = -^—, K'RQ =^ , pp = cos (np, fi).
^PPQ^PRQ
Здесь n(Q,i) — плотность энергии излучения (в единице объема).
В отличие от интегро-дифференциального уравнения (31), (32) уравнения (61), (62) описывают перенос излучения с использованием однородных функциональных средств — интегралов, содержат описание геометрии расчетной области — для пучка пересекающих область лучей, сходящихся в фиксированную точку-цель Q.
3 Ю. А. Дементьев
Глава 3
НОВАЯ ФОРМА ИНТЕГРАЛЬНОГО УРАВНЕНИЯПЕРЕНОСА ИЗЛУЧЕНИЯ
Первая форма интегрального описания распределения лучистой энергии точечного источника
В (41) на основе закона Бугера-Ламберта описан вклад в интенсивность фиксированной точки-цели Q в фиксированный момент времени t — любого из предшествующих точечных источников Р на луче PQ. Объединяя (складывая) вклады всех точечных источников P.R на границе и внутри объема, получаем известное интегральное уравнение Пайерлса, гл. 2.
Новая форма интегрального уравнения переноса тоже получено на основе закона Бугера-Ламберта, но в другой форме (42) — после замены переменной при описании времени. Момент времени t фиксируется для точечного источника Р на луче PQ. Описан вклад этого источника в последующие моменты времени в интенсивность точек-целей Q на этом луче.
Рассмотрим пучок лучей, расходящихся из точки Р и засекающих видимую часть У(Р), S(P) расчетной области, не затененную участками невыпуклой граничной поверхности.
Новое интегральное уравнение переноса описывает распределение энергии, излученной точечным источником Р: распределение проникающей части энергии на граничной поверхности S^P) и распределение выведенной веществом части энергии внутри объема У(Р) вещества, ограниченного поверхностью S.
Распределение энергии описывается:
с учетом видимости — в случае невыпуклости граничной поверхности;
с учетом запаздывания — в случае нестационарности задачи;
с учетом ослабления на основе закона Бугера-Ламберта — вследствие взаимодействия излучения и вещества среды переноса;
почти с теми же весами углового распределения для лучей, что и при записи известного интегрального уравнения (61), (62), но на основе уравнения (42) для луча.
Для граничного точечного источника Р новая форма интегрального уравнения переноса излучения имеет вид:
J+(P,t)dSPdt = | (Jp(Q,t-VpPQ/c)dSQdt-V
+ | ua,sP(R,t + Ppr/c) dVRdt. (63)
В левой части уравнения (63) — исходная величина импульса для граничного точечного источника Р. В правой части — описание распределения этой величины — распределение энергии проникающего излучения в точках видимой из Р части S(P) границы, распределение энергии поглощенного и рассеянного веществом излучения во внутренних точках видимой части V(Р) объема.
Здесь Jp(Q,t + ppg)с) — вклад точечного источника Р в плотность одностороннего потока энергии на границе области, размерность [Jp] = Вт/(м2с); ua,sp(Ri t + Ppr/c) — вклад точечного источника Р в величину плотности выведенной веществом энергии — внутри объема, размерность [na,sp] = Вт/(м3с),. Нижние индексы указывают на принадлежность функции к описанию физических процессов (взаимодействия излучения и вещества) и содержат указание на точечный источник излучения: ua,sp — количество энергии источника Р, поглощаемой и рассеиваемой единицей объема вещества за единицу времени, более подробно — см. во введении. Доли энергии Jp, ua)Sp определяются для лучей PQ, PQ+, PQ — с учетом долей углового распределения, KPQdSQ, KPQ+dSQ+, KpQ-dSQ-, и с учетом ослабления веществом в соответствии с законом Бугера-Ламберта, ехр(—тру), ехр(-тР(3+), ехр(-тр9-):
Jp(Q,t + Ppq/c) dSqdt = J+(P,t) dSpdtKpQdSQexp(-TpQ),
11^a,sP — dVR^a,sP — dVR| Jp(Q+,P)dSQ+- | JpQ-,HdSQ- Ls+s-| Jp(Q+,P)dSQ+- | JpQ-,HdSQ- Ls+s-(64)(64)
1
1
^a,sP — dVR
^a,sP — dVR
| Jp(Q+,P)dSQ+- | JpQ-,HdSQ- Ls+s-
| Jp(Q+,P)dSQ+- | JpQ-,HdSQ- Ls+s-
(64)
(64)
Jp (Q+ ,t+}dSQ+dt = J+(P, t) dSp dtKpQ+dSq+ exp (—tPq+),Jp(Q”,r)dSQ-dt = J+(PP dSpdtKPQ-dSQ- ехр(-тру-),
KPQdSQ = ^dSQ,
^PPQ
KPQ+dSQ+ = ^dSQ+,
^Ppch
KpQ-dSQ- = ^dSQ_.
^PpQ-
Отметим формулу (64) для плотности выведенной энергии излучения веществом элементарного объема dVp. Внутри квадратных скобок — выведенная энергия источника Р внутри объема dVp — разность между энергией, вошедшей в объем вдоль лучей источника Р через освещенную часть 5+ боковой поверхности, и энергией, покинувшей объем через теневую часть S-, — с учетом закона Бугера-Ламберта и весовых функций KPQ+dSQ+, KpQ-dSQ- углового распределения для лучей.
Элементы площади боковой поверхности 5+, 8~ — бесконечно-малые величины, но по ним выполняется интегрирование, рассматриваются элементы их разбиения dSQ+, dSQ-.
При закрепленной сети граничных и внутренних ячеек лучи используются как средство достижения точности при расчете коэффициентов распределения энергии точечного источника, распределения — по ячейкам сети.
Для определения значения подынтегральной функции проводится геометрический расчет пересечения ячеек многими лучами. Это существенная особенность всего рассмотрения. Она возникает при переходе от расчета коэффициентов видимости для переноса в вакуумном приближении к расчету коэффициентов ослабления для переноса с учетом взаимодействия излучения и вещества.
Для определения функции ua,sp(RP + Ppr/с) выведенная элементом объема dVp энергия приводится к единичному объему.
На множитель dt обе части уравнения (63) можно сократить. Однако, удобно рассматривать распределение энергии, а не ее производных.
Для внутреннего точечного источника Р новое интегральное уравнение переноса излучения имеет вид:
uT,s(P, t) dVp dt = | Jp (Q, t + Ppq/c) dSq dt + j n^sP AVr dt;
S(F)(65)
JP{Q,t + Ppq/cPSq = uTiS(J\ t)dVpK'pQdSQ ехр(-трд); (66)
1^a,sP — nz AVr1^a,sP — nz AVrj Jp(Q+,t+)dSQ+ - J Jp^Q",ndSQ- ,s+S-J
1
^a,sP — nz AVr
1
^a,sP — nz AVr
JP(Q+,tpdSQ+ = urtS(P,t')dVpexp(-TPQ+>)KpQ+dSQ+, Jp(Q~ ,t”)dSQ- = ur-s(P .t^dVp exp (-TpQ-^KpQ.dSQ-,
K'PQ+dSQ+ = —dSq+,
KWdS9-= dSQ-■
^PpQ-
Уравнения (63)-(66) дают возможность описания переноса лучистой энергии в соответствии с исходной физической моделью корпускулярного переноса излучения — независимого переноса энергии от каждого точечного источника.
Впервые найдено балансное аналитическое описание функции распределения энергии, излучаемой точечным источником.
Сократим левые и правые части уравнений (63), (65) на величины, описывающие энергию, J+{P,t) dSp dt и иТ,8(РД) dVp dt, соответственно, получим аналитическое описание функции распределения энергии точечного источника: и для проникающей части, и для выведенной части этой энергии:
для граничного точечного источника
1= j exp(-TpQ)KPQdSQ +
S(F)
| exp(_-TpQPKpQ+ dSQ+ -
| exp^-TpQ-^KpQ- dSQ-
s-
(67)(67)для внутреннего точечного источника
(67)
(67)
J exp(-TpQ)K’pQdSQ + S(F)
r J 1
] exp(-TPQPK'pQ+dSQ+ - Ls+] exp(-TPQPK'pQ+dSQ+ - Ls+J \dV^ vp) v
] exp(-TPQPK'pQ+dSQ+ - Ls+
] exp(-TPQPK'pQ+dSQ+ - Ls+
exp^-TpQ-^Kpq. dSQ-exp^-TpQ-^Kpq. dSQ-dVR. (68)dVR. (68)
exp^-TpQ-^Kpq. dSQ-
exp^-TpQ-^Kpq. dSQ-
dVR. (68)
dVR. (68)
Справедливость равенств (67), (68) проверяется непосредственно. Правильность нормировки (л, 4л) угловых коэффициентов обсуждалась во введении. На всех внутренних границах, разделяющих два элементарных объема dVR, dV^ при описании распределения выведенной энергии внутри объема вещества, энергия выходящего потока (через часть теневой поверхности 8~ одного элементарного объема dVR^ равна энергии входящего потока в соседний объем dV^ (через эту же поверхность, но являющуюся теперь частью освещенной боковой поверхности S^ этого соседнего элементарного объема). Но энергия выходящего потока входит в сумму со знаком «—».
После этого замечания в (67), (68) справа после сокращений остаются по два слагаемых. Одно слагаемое (со знаком «+») относится к точке-источнику Р ввода энергии внутрь объема V расчетной области — вершине пучка расходящихся лучей, это слагаемое равно 1. Второе слагаемое (со знаком «—») — интеграл по граничной поверхности S — объединение частей теневых поверхностей S- прилегающих к границе элементарных объемов dVR.
Равенства (67), (68) доказаны:
для расходящегося пучка лучей с центром в точке Р границы S
1= j exp (—tpq^Kpq dSQ + 1 - j ехр^-TpcpKpQdSQ,
S^S^(69)
для расходящегося пучка лучей с центром в точке Р объема V
1= j exp (-tpq)Kpq dSQ + 1 - j exp^-TpQ^KpqdSQ.
s^S(P)z7m
Формулы (69), (70) аналогичны (60).
Сравнивая сказанное относительно (67), (68) и ранее — относительно (54), можно говорить об аналоге формулы Ньютона- Лейбница, (69), (70). Формула (54) описывает перенос и ослабление вдоль одного луча. Формулы (67), (68) описывают перенос и ослабление вдоль пучка лучей с центром в точечном источнике Р излучения, граничном или внутреннем, соответственно. Подобно формуле Ньютона-Лейбница интегралы, описывающие
распределение на луче выведенной энергии после сокращения величин, относящихся к внутренним поверхностям объема, сводятся к двум членам суммы, первому и последнему.
Вторая форма интегрального описанияраспределения лучистой энергии точечного источника
Уравнения (63)-(68) сложнее численной схемы.
Рассмотрим второй вариант интегрального описания переноса излучения точечного источника в соответствии с фотометрической моделью — аналоги формул (63), (65), (67), (68) на другой аналитической основе: на основе интеграла вдоль луча — на основе формулы интегрирования Ньютона-Лейбница и интеграла Стилтьеса. Опираясь на (54) для функции распределения вдоль одного луча выведенной веществом части излучения точечного источника и формулу (55) описания распределения излучения точечного источника вдоль луча, получаем:
для расходящегося пучка лучей с центром в точке Р границы S
Q
dE^P, t) = dE(P, t) | ехр(-трд+^-ехр(-тря-) + S(P) p
+ exp (-tpq) dSq =
zQ.
= dE(P,t) j H d(-exp [-r(g)]) + exp(-rpQ) \KPQdSQ,
S(P) Ц
(71)
^a^PQ^R' ^ + PPR) AVr dt =
= dE(P,t)KPQdSQ\exp(-TpR+) - exp (-трд-)],
Jrq(Q, t + ppg) dSQ dt = dE(P, t)KPQdSQ exp (-tPq),
KPQdSQ = (mpMq)/(^pq)^5'q-
для расходящегося пучка лучей с центром в точке Р объема V
Q
dE(P, t) = dE(P, t) | ехр(-тдд+)^-ехр(-тдд-) + S(P) P
+ exp pTpq^ }kPq dSq =
= dE(P,t) j j d(-exp[-r(g)]) + ехр(-тР$) VKpqdSq, S(P) ^P
(72)
■U-a.sPQ^R, t + Ppr) dVp dt =
= dE(J\t)lPpQdSQ[exp^^ - exp (-тРД-)],
Jpq(Q, t + Ppq) dSq dt = dE^P, t)K'PqdSq exp (-tPq),
KpqdSq = pqf^ixppq^dSq.
Рассматриваются расходящиеся лучи PQ, ведущие из точки- источника F, граничной или внутренней, в точки-цели Q границы S с учетом видимости, Q е S^P).
Здесь dE(P,t) — импульс, излучаемый за время dt точечным источником Р, граничным или внутренним, соответственно, — за счет плотности одностороннего потока окрестности границы, dE(P,t) = J+{P,t) dSp dt, или за счет излучения вещества окрестности объема, собственного излучения и рассеяния, dE(P, t) = uT,s(P, t) dVp dt. Мы решаем нестационарную краевую задачу, нам необходимо локальное описание точечного источника, с локальным определением точки и окрестности границы или вещества, с локальным указанием момента времени и диапазона времени излучения. Сравните сказанное с функциями описания источника правой частью интегро-дифференциального уравнения.
dE(P,t)KpQdSQ или dE(P,t)KpgdSQ — результат углового распределения энергии, часть импульса, излучаемая в направлении PQ внутрь телесного угла dwpQ или du'pg, — внутрь конической поверхности с вершиной Р (граничного или внутреннего источника), опирающейся на окрестность dSq точки-цели Q границы S.
Мы исчерпали функции правой части интегро-дифференциального уравнения, описали источник излучения и его интен-
сивность (с множителями) в выбранном направлении — тоже локально, локально по телесному углу, не на единицу телесного угла, как в (38),
dE(P, t)KPQdSQ = Р(Р, t)dSPdtnp^2-, (73) j+ = тгЦРД) — закон Ламберта,
dE(P, t)K'PQdSQ = ur.s(P, t^dVpdt1^., (74)
Ur,8 = Xa,ScU = И^РРхрРр, (16), (18).
Для описания углового распределения излучения точечного источника — на единицу площади граничной поверхности, нами использованы Kpq, KpQ, и тем самым уточняются возможности интегрального уравнения Пайерлса. В нашем случае описывается локальное распределение, KpqdSQ, KpqdSQ, — внутрь конуса dupQ или dupQ, а не в сторону единичной площади граничной поверхности.
Далее в (71), (72) используется решение характеристического уравнения, закон Бугера-Ламберта в форме (42):
R — центральная точка элемента объема dVp, пересекаемого лучом PQ; dE(P,t)KpQdSQexp[—r(PR+)] или dE^P^KpqdSq exp [—т(РР^У\ — энергия проникающего излучения, достигшая вдоль луча PQ освещенной стороны боковой поверхности объема dVp; dE^P^KpqdSq ехр[—т(РД-)] или dE^P^KpqdSq exp [—t(PR-^ — энергия проникающего излучения, достигшая вдоль луча PQ теневой стороны боковой поверхности объема dVp;
разность — энергия, выведенная веществом объема dVp из конуса dcopq на отрезке пересечения dgp = R4 R.
В (64) и (66) выведенная элементом объема dVp энергия — для излучения фиксированного точечного источника, вычислялась целиком и делилась на dVp.
Необходимость деления на dgp в (71), (72) подтверждается проверкой размерности. Интегрирование вдоль луча по формуле Ньютона-Лейбница как раз то, что нам нужно. Эта удивительная формула означает сложение разностей без умножения каждой из них на элемент интегрирующей величины. Сокращение средних слагаемых характерно и для формулы Ньютона- Лейбница, и для распределения выведенной энергии вдоль луча.
Все сказанное означает необходимость интегрирования характеристического уравнения. Из уравнения следует: производ-
74 Гл. 3. Новая форма интегрального уравнения переноса излучения ная от интенсивности равна линейной плотности выведенной энергии.
Будем называть Ua^pQ^R^t) лучевой долей распределения плотности выведенной энергии, размерность [i^pq] = Вт/(м3с);
dE(P,t)KpQdSQexp(—TpQ) — часть импульса dE(P,t), достигшая границы вдоль луча PQ; размерность [dE(P,t)KpQ х х dSq exp (-tPq)] = Вт;
будем называть Jpg(Q, t + Ppq) лучевой долей распределения плотности энергии проникающего на границу излучения, размерность Upq(Q^ + PPQ^ = Вт/(м2с).
Вторая форма интегрального уравнения целиком соответствует численной схеме.
Так как конусы dupQ, du/pg заполняют (см. введение) всю видимую часть объема расчетной области, сумму лучевых долей плотности выведенной энергии ua,sPQ(R,t) с фиксированным R назовем вкладом ua,sp(R, ^) точечного источника Р в величину плотности выведенной энергии ua)S(R,t\
Конусы dcupQ, du/pg всех точечных источников многократно заполняют всю расчетную область и его часть — элемент объема dVp. Вклады ua,sP(R, i) всех точечных источников, граничных и внутренних, одновременно приходящие в центр R объема dVp, складываются в величину плотности выведенной энергии иа^РД) — веществом объема dVp.
Так все и будет в численной схеме.
Излучаемая импульсами энергия, (73), (74), распределяется в шаровых слоях толщины dg = cdt (по радиусу) с центром в точке Р источника; в слое (0, dg) распределена энергия импульса dE(P,t)KfRgdSQ, излученного в текущий момент времени; в слое (5dg, 6dg) распределена энергия импульса dE{P,t — 5dt)KfRQdSQ, излученного ранее — на пять шагов dt по времени.
Проверка правильности формул (71), (72) выполняется на основе формулы интегрирования Ньютона-Лейбница — сокращением средних слагаемых.
Деля все члены уравнений (71), (72) на общий множитель dE(P,t), находящийся за пределами интегрирования, получаем описание функции распределения единичного импульса, или функции распределения вероятности для излучения точечного источника:
для расходящегося пучка лучей с центром в точке Р границы S
1 = j p(-exp[-T(9)])+exp(-rpQ) KPQdSQ, (75) S(F) ^FJ
для расходящегося пучка лучей с центром в точке Р объема V
= | Н dp exp [-т(у)]) + exp (-tpq) >KpQ dSQ.
S(P)1
Уравнения (71)-(75) значительно проще уравнений Пайерлса (61), (62) и первой формы уравнений (64), (66). В уравнениях (64)-(66) и в уравнениях Пайерлса в правой части уравнений содержатся по два интегральных члена — для описания переноса на границу и для описания переноса внутри объема.
Уравнения (71)-(75) описывают перенос для пучка расходящихся лучей, аналитически описывается распределение энергии вдоль фиксированного луча PQ от точечного источника до граничной точки в области видимости источника: ввод энергии точечного источника в конус телесного угла, пересечение лучом объема поглощающего и рассеивающего вещества, выход проникающей части излучения на границу.
3.3. Обсуждение результатов
Коснемся особенности описания углового распределения излучения точечных источников в (63)-(75), в отличие от аналогичного описания в уравнении Пайерлса.
Множители KpqcISq, KpgdSq в (71)-(75) содержат в себе описание телесных углов dujpq, du/pg — конусов с вершиной в точечном источнике Р, граничном или внутреннем, опирающихся на элемент граничной поверхности dSg,
KpQdSQ = dE^EQ, K'pQdSQ = ^.(76)
В уравнении Пайерлса (61), (62) множители Kpq, KpQ, определяющие угловое распределение излучения в точках-источниках, не включают в себя элементарные величины площади поверхности в окрестности точки-цели, как dSg в (76). Элемент поверхности dSp относится в (61), (62) к описанию источника излучения. Угловое распределение в (61), (62) связывает источник
излучения с единицей поверхности в окрестности точки-цели. Для более детального описания углового распределения в уравнениях Пайерлса прибегают к повторному интегрированию по поверхности S, чтобы ввести в рассмотрение dSq. Угловые коэффициенты в методе счета коэффициентов видимости определяются многократными интегралами на основе уравнения Пайерлса, например в [17]. При использовании уравнений (71)-(75) все коэффициенты определяются лишь двукратными интегралами — интегралами по телесному углу.
Переход от интегралов к любой допредельной частичной дискретной сумме делает сами уравнения (71)-(75) непосредственно реализуемой численной схемой. Промежуток dt времени становится шагом по времени.
На основе фиксированных координат точек Р, Q вычисляется коэффициент углового распределения KpqdSq или KpgdSq для луча PQ. После умножения углового коэффициента на энергию импульса J+{P,t) dSp dt или uT)S(JM) dVp dt определяется стартовая энергия, высвечиваемая точечным источником внутрь телесного угла вокруг луча PQ.
Схема численного описания переноса с помощью фотометрической модели на основе уравнений (71)-(75) называется согласованной схемой. Все доли распределения импульса энергии точечного источника вычисляются согласованно с отслеживанием последовательности углового распределения и ослабления энергии при движении излучения вдоль луча. Сумма всех долей распределения строго равна распределяемому импульсу. Уравнения (75) используются для вычисления коэффициентов распределения; формулы (71), (72) в дискретном варианте используются для описания переноса (22).
Ер^ = У^ Ер^^т + У^ Ep^ij, У^ Ргт 4" У^ Pij = 1 •m 3 т з
Лучевые доли распределения плотности энергии Jpg^Q^t) проникающего в момент времени t в граничную точку Q с окрестностью dSq излучения, (71), (72), одновременно приходящие от разных точечных источников Р, граничных и внутренних, складываясь, образуют полные величины J-(Q,t) плотности одностороннего потока проникающего излучения,
ГШЖг = j JpQ(Q,t)6QpdSP+ | JpQ(Q,tWQpdVP,S(Q) V(Q)
где dqpdSp = dSq, dqpdVp = dSq.
Лучевые доли распределения плотности ua,sPQ(R,t) выведенной в момент времени t веществом внутренней точки R с окрестностью dVp энергии, (71), (72), одновременно приходящие от разных точечных источников Р, граничных и внутренних, складываясь, образуют полные величины ua,s(R, i) плотности выведенной энергии,
Ra,s{RP}dVp = | | Uu spq^R, t^dppQ dSp dSq +
S(R) S(R)
+ j | Ua>sPQ<R,t^^pQdVpdSQ, V(R) S(R)
где 5^pQdSpdSQ = dVR, 6^PQdVpdSQ = dVR.
В дискретном случае накопление величин суммированием лучевых долей от разных точечных источников выполняется с учетом запаздывания со вторым порядком точности по времени благодаря специальному алгоритму размещения лучевых долей в матрицах переноса ещ ew (гл. 4) — с определением целей переноса и групп запаздывания.
Можно было бы вычислять выведенную веществом энергию для каждой дискретной ячейки объема отдельно, следуя непосредственно формуле (64) — находить вклад в эту энергию для каждого источника по каждому дискретному лучу. При этом для определения оптической толщины пути до входа, необходимо повторять вычисления координат точек пересечения лучей с границами других предшествовавших элементарных объемов dVp.
Назовем базовым множеством лучей множество лучей PQ, ведущих в счетные центры Q элементов dSq разбиения граничной поверхности. Базовое множество лучей определяет базовое разбиение телесного угла для точки-источника Р.
Согласованным счетом всех интегралов по телесному углу (2л) для граничных источников или (4л) для внутренних источников будем называть их вычисление с использованием только базового множества лучей и базового разбиения телесного угла. При этом используется «наведенное» разбиение боковых поверхностей 5+, 8~ пересекаемых элементов объема dVp с вычислением точек пересечения R"^ R , см. рис. 1.
По каждому лучу базового множества вычисления проводятся один раз: от точки-источника Р до точки Q на границе области. Вычисляются координаты точек входа R^ и выхода R — для пересекаемых лучом элементарных объемов dVp. Определи-
ется энергия, выведенная веществом этого объема. И т. д., нако
нец оставшаяся энергия становится слагаемым при вычислении граничной величины Jp(Q,t + Ppq/с).
При таком согласованном выборе дискретных средств для
приближенного вычисления определенных интегралов должно
Согласованный расчет коэффициентовСогласованный расчет коэффициентовРис. 3.1. Согласованный счёт коэффициентов Здесь Р — точечный источник, Q Е S, R Е с/цд; 8"^, S- — видимая и невидимая стороны боковой поверхности элементарного объёма dvR для точки Р;Q^ Е 8^, Q- Е 8"Рис. 3.1. Согласованный счёт коэффициентов Здесь Р — точечный источник, Q Е S, R Е с/цд; 8"^, S- — видимая и невидимая стороны боковой поверхности элементарного объёма dvR для точки Р;Q^ Е 8^, Q- Е 8"ванного луча. Используютсяванного луча. Используютсябыть выполнено требование стремления к нулю диаметров областей разбиения. Оно будет выполнено, если разбиение граничной поверхности и внутренних дискретных объемов проводится — даже для двумерных расчетов — по всем трем координатам. Например, для двумерных расчетов — по цилиндрическим координатам, включая угол вращения.
Согласованный расчет коэффициентов
Согласованный расчет коэффициентов
Рис. 3.1. Согласованный счёт коэффициентов Здесь Р — точечный источник, Q Е S, R Е с/цд; 8"^, S- — видимая и невидимая стороны боковой поверхности элементарного объёма dvR для точки Р;
Q^ Е 8^, Q- Е 8"
Рис. 3.1. Согласованный счёт коэффициентов Здесь Р — точечный источник, Q Е S, R Е с/цд; 8"^, S- — видимая и невидимая стороны боковой поверхности элементарного объёма dvR для точки Р;
Q^ Е 8^, Q- Е 8"
ванного луча. Используются
ванного луча. Используются
Выше и во введении мы говорили о квазианалитичности и контролируемой точности описания углового распределения импульса энергии. Отметим аналитическое описание (71)-(75) распределения излучения точечного источника вдоль фиксиро- решение (42) характеристическо
го уравнения и формула интегрирования Ньютона-Лейбница. На этой основе совершенно точно — со всеми вычисляемыми ЭВМ знаками — определяются доли распределения выведенной веществом энергии (разность экспонент) и доля проникающего на границу излучения. Для сравнения и в качестве упражнения
предлагаем читателю подумать над аналитическим определением (например, с помощью интегрирования) лучевой доли плотности излучения Urq^RR + ррр/с) и вычислением среднего расстояния ррр. Нужно учитывать особенности поведения экспоненты, в данном случае, — на отрезке dgp = R+ R~. Экспонента отличается от линейной функции большей степенью роста (в данном случае — в сторону точки R^), чем любая степенная функция аргумента д. Более быстрым ростом обладают экспоненты с большей степенью и далее — экспоненты со степенью — экспонентой, типа ехр [—ехр (р)]. Плотность излучения U(RR + ppp)^ не входит в число опорных функций, определяющих точность вычислений в нашей схеме описания переноса. Опорными функция-
ми, вычисление которых необходимо для вычисления следующих функций, на шаге по времени являются: — для граничного точечного источника F-J+(F,£), ua,sPQ(R, tPR), J-^QPPqY — для внутреннего точечного источника Р — ur,s(P,t), ua,sPQ(R,tPP), иа)8(РДРР). В других методиках U(R,t + рРя/с) — опорная функция.
При решении краевой задачи для замыкания вычислений граничных величин используются граничные условия (9)-(11), (94) [7].
Для замыкания вычислений величин, описывающих перенос внутри расчетной области, используется уравнение локального баланса энергии для вещества элемента объема dVPl [8]:
dp •
= Ua-Ur.(77)
dt
Сначала на основе принципа Шустера выведенная энергия uaiS разделяется на поглощенную и рассеянную энергию, ^a,s^a^a,s^s
Ua = —, Us = .(78)
При решении нелинейных уравнений (20), (94) относительно неизвестной величины — температуры вещества — используются локальные итерации Ньютона, как и в других методиках приближенного описания переноса излучения, [8].
Вблизи граничной поверхности имеет место разрыв решения в связи с двумя, вообще говоря, разными значениями температуры источников излучения:
эффективной температурой вещества за пределами расчетной области, найденной при решении уравнения (94) граничного условия,
температурой вещества изнутри расчетной области, найденной при решении уравнения (20), (77).
В гл. 1 мы уже касались описания этого разрыва в численных схемах на основе интегро-дифференциального уравнения. Многие эффекты переноса описываются на уровне построения схемы счета, минуя прямую аппроксимацию уравнения переноса.
На самой граничной поверхности могут быть разрывы температуры [12], тогда из такой точки происходит перенос трех разнотемпературных интенсивностей. Разрывы решения возникают и в точках разрыва коэффициентов взаимодействия или на границах заданной функции источника. В численных методиках такие эффекты учитываются при выборе счетной сети ячеек.
Большую роль в развитии теории описания излучения нагретого тела сыграло рассмотрение Г. Р. Кирхгофом равновесно-
80 Гл. 3. Новая форма интегрального уравнения переноса излучения го излучения. Внутри области с помощью граничных условий и начальных данных установлен режим изотермии, температура T(F,t) постоянна для всей области и для границы. Для уравнения теплопроводности, например, в таких условиях правая часть уравнения равна нулю — отсутствует градиент температуры — и, в результате, температура постоянна и по времени.
А вот модель переноса излучения должна приводить к тому же результату совсем другим образом. Рассеяние в этом расчете не нужно рассматривать. Стационарный режим должен установиться как режим равновесия — между процессами поглощения и высвечивания. Точечные источники высвечивают импульс энергии, соответствующий температуре. Из-за ошибок в угловом распределении вдоль каких-то лучей пойдет завышенный поток энергии, вдоль других лучей поток энергии будет сбалансированно занижен. Расчет состояния изотермии, по-видимому, — серьезный тест для методик численного описания переноса излучения. Чем больше неточностей в деталях, тем больше отклонений от постоянной температуры. По данной методике такой расчет не проводился.
Глава 4
ЧИСЛЕННАЯ МОДЕЛЬНЕСТАЦИОНАРНОГО ПЕРЕНОСА ИЗЛУЧЕНИЯВ ДВУМЕРНОМ ПРИБЛИЖЕНИИДЛЯ ШИРОКОГО КЛАССА СЕТОК
Сеть шагов по времени.
Сеть граничных и внутренних ячеек
Счет нестационарных краевых задач будем вести шагами по времени к преимущественно с постоянным шагом т по времени.
Введем сеть граничных и внутренних ячеек.
Граничные ячейки и их площадь будем обозначать S^ S3. Нижние индексы — номера граничных ячеек, 1 < i,j < kgt, kgt — число граничных ячеек.
Внутренние ячейки и их объем будем обозначать Vmi Vn. Нижние индексы — номера внутренних ячеек, 1 < m,n < k^t, kvt — число внутренних ячеек.
Первые индексы (г или т) будут, как правило, использоваться для указания номера ячейки — источника излучения, вторые (у или п) — для указания номера ячейки — цели переноса.
В случае решения двумерных задач предполагается, что расчетная область имеет ось симметрии, граничные ячейки — поверхности вращения, внутренние ячейки — объемы вращения.
Разбиение расчетной области на ячейки в этом случае не удовлетворяет требованию к разбиению области интегрирования при определении интеграла как предела частичных сумм: диаметр элементов разбиения — поверхностей и объемов вращения — не стремится к нулю при увеличении числа разбиений.
Поэтому дополнительно разобьем сеточные ячейки-цели, граничные и внутренние, равномерно по углу вращения. Обозначим полученные доли разбиения поверхностей и объемов вращения: Sjki Vnb 1 < к < к^ где fc^ — заданное число равномерных разбиений по углу вращения ip.
Фактически при решении двумерных нестационарных задач мы прибегаем к трехмерной сети ячеек-целей, в частности, для
правильного описания эффектов запаздывания. Расстояния от фиксированной ячейки-источника до долей разбиения ячейки-цели разные. Двумерными будут функции, описывающие результат переноса. За время фиксированного шага по времени в ячейку- цель будет перенесена энергия, равная сумме долей энергии, перенесенных от всех источников во все доли разбиения этой ячейки по углу — в течение этого времени.
С помощью выбора сети граничных и внутренних ячеек описывается состояние вещества и распределение источников излучения. Сеточные точечные источники расположены в счетных центрах Р^ Рт ячеек, граничных и внутренних — в полуплоскости % = О, Pi G Si, Рт
Основные счетные величины описания веществаи излучения на шаге по времени
Основными счетными величинами на границе будут величины плотности односторонних потоков лучистой энергии J^", Jr — высвечиваемой внутрь расчетной области и приходящей на границу со стороны расчетной области, соответственно, [2-6]. Верхний индекс (+) соответствует направлению вектора внутренней нормали к границе по отношению к расчетной области. Индекс (—) — обратному направлению — в окружающее пространство.
Размерность плотности односторонних потоков \J^, Jr] = = Вт/(м2с).
Односторонний поток J^ определяет, (8), импульс ET^t) энергии излучения, высвечиваемой ячейкой Si за шаг по времени из граничной точки-источника Pi внутрь расчетной области, E^(t) = J^(t)SiT. При описании переноса импульс ET^t) распределяется по направлениям счетных лучей PiQj неизотропно — пропорционально косинусу угла наклона вектора источник-цель PiQj к вектору внутренней нормали пг к граничной поверхности S (закон Ламберта) — и пропорционально доле телесного угла Awij, отнесенной к этому лучу.
Односторонний поток Jr определяется в ходе расчета переноса, импульс Ej(ty Ej(t) = Jr(t)SjT, объединяет проникающую часть энергии от всех точек-источников, граничных Pi и внутренних Рт.
Предполагается, что односторонние потоки J^, Jr связаны граничным условием (9), J^" = f(JjQ, позволяющим замкнуть вычисления в граничных ячейках — рассчитать энергию, излуча-
емую граничными источниками, J^ на новом шаге по времени — на основе энергии, J^ перенесенной на границу на предыдущем шаге по времени.
Основными счетными величинами во внутренних ячейках описания вещества будут величины плотности энергии иг,8т = — Wm ~I- ^snu ^a,sn — ^ап ~I- ^sn^ соответственно, излученной внутрь расчетной области или выведенной веществом, — в единицу времени. Нижние индексы поясняют:
г — за счет собственного излучения вещества,
s — за счет рассеяния излучения веществом, а — за счет поглощения.
Размерность плотности излученной, выведенной, поглощенной или рассеянной энергии \uT)Smi ua,sn, иап, uSTl\ = Вт/(м3с).
Величина итт определяет лучеиспускательную способность вещества точки-источника, описывается энергия собственного излучения вещества — пропорционально четвертой степени температуры вещества — в условиях локального равновесия излучения и поглощения вещества, это законы Кирхгофа и Стефана- Больцмана.
Величина и8т определяет способность вещества к рассеянию излучения, описывается энергия квантов излучения, изменивших направление переноса внутри данной ячейки при ее пересечении лучами, идущими от всех сеточных источников, — на предыдущем шаге по времени.
Величина ur,sm определяет возможности внутреннего сеточного источника излучать энергию на шаге по времени — в соответствии с температурой вещества — и рассеивать излучение — по результатам счета в этой ячейке на предыдущем шаге по времени.
При описании переноса излучения предполагается, что энергия иг,8т собственного излучения вещества и рассеянная энергия распределяются по направлениям лучей изотропно, т. е. независимо от направления луча. Энергиявнутренних точечных источников распределяется по лучам PiQj пропорционально лишь долям телесного угла А^, отнесенным к этим лучам.
Величина иап определяет способность вещества пересекаемой лучом ячейки к поглощению излучения.
Величина и8П определяет способность вещества пересекаемой лучом ячейки к рассеянию излучения.
Величина иа)8Т1 определяет способность вещества пересекаемой лучом ячейки к поглощению и рассеянию излучения —
84 Гл. 4. Численная модель нестационарного переноса излучения способность вещества к выводу энергии из пронизывающего его излучения.
Каждой сеточной ячейке ставится в соответствие точечный источник излучения.
Предположим, что на предыдущем шаге нами вычислены во всех граничных и внутренних ячейках, соответственно, плотности энергии J^t) и uT)Sm(tY излученной в единицу времени внутрь расчетной области. Тогда в текущий момент времени t определим, по (7), (8), величину импульса энергии, EJ(t) или E^ty излученного каждым сеточным источником внутрь расчетной области за шаг т по времени:
Е^Ц = Г(Ц8гТ, Е^Ц = uT,sm(t)VmT. (79)
Нижние индексы в (79) поясняют физический смысл величин и номера счетных источников и целей.
В предположении квазистационарности геометрии расчетной области и величины пробегов излучения квазианалитически вычисляются коэффициенты распределения излучаемой энергии. Точечному сеточному источнику ставится в соответствие дискретный пучок расходящихся счетных лучей. Для каждого луча вычисляются по очереди точки входа и выхода из пересекаемых ячеек, оптическая толщина пути и ослабление излучения на основе закона Бугера-Ламберта. Нестационарный перенос энергии описывается со вторым порядком точности учета запаздывания по времени.
Отложим пока обсуждение алгоритма вычисления коэффициентов и алгоритма учета запаздывания по времени.
Распределение импульса излучения сеточногоисточника
В результате переноса проникающая часть энергии Е^ или Е^, излученной за шаг т по времени каждым сеточным источником, дойдет до границы и распределится по ней, другая — выведенная веществом — часть энергии распределится внутри объема вещества:
(^) = У? -^ijk^ + Pijk/+ У^ -^гпк^ + Pink/с\ (80) jkпк
^mto = У? ^mjk^ + Pmjkf + У^ Emnk(t + pmnk/cY (81) jkпк
3Десь Е^к, Е{пк - доли энергии Е^ E^jk, Е^пк - доли энергии Е^.
Обозначим эти доли:
■^ijk = Pijk^i , -^ink = Pink^i i
Т _ т т _ т-^mjk PmjkEm1 -^mnk PmnkEm-
Обозначенные доли назовем базовыми коэффициентами ослабления. При описании переноса будут использоваться коэффициенты переноса, вычисленные на основе базовых коэффициентов ослабления.
Внутри скобок в (80), (81) указаны: момент времени t излучения импульса энергии — момент времени текущего шага — и моменты времени достижения цели долями импульса — с учетом запаздывания. Время запаздывания определяется средними расстояниями р^у, р^у, pmjk, между ячейкой-источником излучения, Si или Vmi и долями (с номером fc) разбиения по углу вращения ячейки-цели, Sj или Vn. Средние расстояния также входят в число базовых коэффициентов. Без дополнительного разбиения по углу вращения поверхностей вращения Sj и объемов вращения Vn разбиение границы и объема расчетной области на ячейки не удовлетворяет требованию стягивания к нулю диаметра областей разбиения при численном интегрировании.
Уравнения (80), (81) являются аппроксимацией уравнений (75).
Все коэффициенты вычисляются до счета переноса. Счет шага по времени разделяется на два этапа: счет коэффициентов и счет переноса.
К вычислению коэффициентов обратимся позже. Будет обсуждаться счет базовых коэффициентов, (см. рис. 2), и счет на их основе коэффициентов переноса с учетом запаздывания. Отвлечемся временно и от учета запаздывания.
Доли Е^к, Е^к проникающей части лучистой энергии, одновременно приходящие на границу, складываясь, позволяют вычислить плотность Jy потока лучистой энергии, приходящего в граничную ячейку Sj со стороны расчетной области.
Точно так же доли Е^., E^nk выведенной веществом части лучистой энергии, одновременно приходящие во внутреннюю ячейку, складываясь, позволяют вычислить плотность ua,sn выведенной веществом внутренней ячейки Vn.
Граничные ячейкиВнутренние ячейки
(в направлении переноса): S„ s,к.
Направление переноса Базовые коэффициенты a? jРу. еу, р,
^?Pm?^m?
Pm? Pm
Рис. 4.1. Схема связи счётных направлений и базовых коэффициентовослабления для излучения граничных источников
Замыкающие вычисления
Замыкающие вычисления — локально в каждой счетной ячейке — аналогичны замыкающим вычислениям [8].
В граничных ячейках на основе граничных условий (9)—(11), (94) (см. ниже) по найденной величине одностороннего потока JC" вычисляется величина другого одностороннего потока J^"; например, с помощью заданной величины aj альбедо,
j+ =.
В каждой внутренней ячейке замыкающие вычисления ведутся — на основе аппроксимации уравнения локального баланса энергии вещества [8], (20),
^ —— ^rni(33)
поглощенная энергия иа определяет приращение за шаг по времени внутренней энергии е^ вещества и величину энергии ит собственного излучения вещества.
Сначала выведенная веществом этой ячейки энергия ua,sn разделяется на поглощенную иатг и рассеянную энергию и8тг.
Используется принцип Шустера, (78),
^а^п^а/п^a^n^sn
Нап = , ^sn = •(84)
^a,sn^a,sn
После вычисления плотности иаи поглощенной веществом энергии итерациями Ньютона [8] решается нелинейное уравнение баланса энергии вещества, (20); (83), — определяется новое значение температуры вещества:
СупРпТп — ^гпп = ^ап ^тп^Т, Urn — ЖапСаТп.
Напомним, что величины энергии uaiSUl иап, usn, uTiSm относятся к единице объема и времени.
После выполнения пересчета по формулам (35)-(37) определяем величины совокупной плотности энергии ur,sm, излучаемой источниками на новом шаге по времени,
^r,sm — ^тт 4" U8m.
С помощью формул (83), (84) счет переноса на одном шаге по времени завершен, мы готовы к повторению счета на новом шаге по времени.
Отметим: все формулы почти одинаковы для счета двумерных и трехмерных задач — для квазидвумерных сеток в локальных дополнительных системах координат [18]. Программы счета коэффициентов и переноса лучистой энергии по методу коэффициентов ослабления для двумерных и трехмерных задач, по-видимому, можно объединить.
Алгоритм учета запаздывания
Формулы (79) для каждого адреса источника определяют дискретизацию излученной энергии импульсами в соответствии с шагами по времени. Формулы (79) описывают сеточное распределение импульсов энергии излучения с указанием адресов источника и цели. Если бы расстояния р^ь ptnb Pmjk, Ртпк между сеточными источниками и целями были бы кратными длине пути ст света за шаг по времени, то и приход лучистой энергии в ячейки-цели был бы полностью согласован с шагами по времени. Диаграмма, представленная на рис. 3, поясняет особенности согласования прихода энергии в ячейки-цели в общем случае, когда расстояние источник-цель не кратно ст,
Pijk = CT(riyjy + C^ijk^y
Рис. 4.2. Диаграмма обмена счётными квантами для граничных ячеек. Pijk = ст(п^к + otijkY в данном случае тг^к = 1
Мы ввели обозначения для целой и дробной частей деления расстояний на ст, соответственно: п^у и
Прямоугольниками с длиной основания ст на рис. 3 изображены два последовательных счетных импульса, ЕГк, Е^к, высвеченной энергии граничным источником Si в сторону элемента (с номером к) разбиения граничной цели Sj.
За текущий шаг по времени излучаемая энергия ЕГк с постоянной интенсивностью заполняет в пространстве слой толщиной ст по нормали к поверхности источника Si на пути к поверхности ячейки-цели Sjk- Для диаграммы выбран случай ^ijk = 1» Pijk = Ст(1 + CXijy^.
Часть (1 — а^к) энергии импульса Е'-к, излученного на предыдущем шаге, к концу текущего шага поглощена поверхностью-целью Sjk- Кроме того, за текущий шаг поглощена часть энергии импульса Е^к, высвеченного на еще более раннем шаге. Этот импульс не изображен прямоугольником на диаграмме.
На следующем шаге будет поглощена часть (1 — а^у) энергии импульса ЕГк и часть а^у энергии импульса ЕЕ.
В общем случае энергия импульса ЕГк будет поглощена частями на двух шагах по времени.
Часть (1 — а^к) энергии импульса будет поглощена через п^у шагов по времени.
Другая часть а^у — через (п^- + 1) шаг по времени.
Аналогично для импульсов EJnk, E^jk, Е^пк.
Подставим обозначения (82) для долей распределения импульса в (80), (81) и сократим обе части полученных равенств на величину импульса сеточного источника ЕТ или Е^. Получим
дискретный аналог равенств (75)
U(^) = У^ PijkA + pijkfA + У^ Ргпк^ + Ргпк/с), (85) jkтгк
Im(^) = У^ PmjkA + pmjk/А + У^ PmnkQt + Ртпк/с),(86)
jkпк
определяющие сеточное распределение импульса источника, соответственно граничного или внутреннего, независимо от величины энергии.
Мы получили описание распределения проникающей части лучистой энергии по граничной поверхности и распределения выведенной веществом части лучистой энергии внутри объема расчетной области.
Вместе с тем, (85), (86) — распределение по времени прихода для долей импульса, излученного в момент времени t.
Базовые коэффициенты и коэффициентыпереноса
Введем понятие коэффициентов переноса в отличие от базовых коэффициентов ослабления: рСк = (1 - ощкУРць п^к, р"-к =
^ijkPijk-
Аналогично: P^, nink, p'kc Pmjk’ nmjk, Pmjk’ Ртпк’ Птпк, Ртпк"
Введенные коэффициенты рассчитываются с помощью базовых коэффициентов и позволяют описать распределение импульсов согласованно с шагами по времени с учетом запаздывания в общем случае, когда вычисленные средние расстояния не кратны ст — со вторым порядком точности по времени:
W) = ^Pijk^t + nijk^ + р'^ + TkijkT + т)} + jk
+ ТМпк^ + ПгпкТ) + р"пкЦ + TlinkT + T), тг к
lm(<) = ^VPmjk^ + nmjkT) + P^jk^ + nmjkT + r)} + jk
+ ^APmnkA + Птп^) + PmnkA + ^mnkT + тгк
Переносимые счетные импульсы, согласованные с шагами по времени, описываются так:
Ej,jk PijkEi iEU=^mkELE,m3k=PmjkELEmnk PmnkEmiEj,jk PijkEi iEU=^mkELE,m3k=PmjkELEmnk PmnkEmiEjjk PijkEi ’Eink PinkEi 'Етпк РтпкЕтлEjjk PijkEi ’Eink PinkEi 'Етпк РтпкЕтл
Ej,jk PijkEi iEU=^mkELE,m3k=PmjkELEmnk PmnkEmi
Ej,jk PijkEi iEU=^mkELE,m3k=PmjkELEmnk PmnkEmi
Ejjk PijkEi ’Eink PinkEi '
Етпк РтпкЕтл
Ejjk PijkEi ’Eink PinkEi '
Етпк РтпкЕтл
Чтобы «доставить» до цели в нужные моменты времени вычисленные доли проникающей части энергии Е^к, Е^к, Efmjk, ЕгЕк, используется массив ev [18]; чтобы «доставить» до цели вычисленные доли выведенной веществом части энергии Efink, Е^, Efmnk, Е^пк, используется массив ew. Линейные массивы ещ еии имеют структуру двумерных массивов. Переносимая энергия в этих массивах распределена по группам запаздывания, а внутри группы — по целям. Группу запаздывания в массиве ег) составляют суммарные импульсы энергии, переносимые одновременно (в течение одного шага по времени) к граничным ячейкам Sj. Группу запаздывания в массиве вии составляют суммарные импульсы энергии, переносимые к внутренним ячейкам Vm в течение одного шага по времени.
Импульсы энергии, составляющие первую группу запаздывания массива ev (или еггггУ достигнут целей Sj (или ¥т) на следующем шаге по времени. Кванты энергии, составляющие вторую группу, достигнут целей через два шага по времени, и т. д.
Импульсы энергии Е’-к будут объединены суммированием с 7-ми числами п^у группы запаздывания массива ег). Импульсы энергии Е^.к будут объединены суммированием с у-ми числами {nijk + 1) группы запаздывания массива ег). И т. д.
На следующем шаге по времени величины энергии, составляющие первую группу запаздывания массива ещ перейдут в массив JJ величин плотности одностороннего потока в новый момент времени, после деления на ^зТУ
Величины энергии, составляющие первую группу запаздывания массива еии перейдут в массив величин плотности ua^sn, после деления на (У^тУ и станут основой для вычисления энергии uTiSTl излучения для внутренних источников.
Оставшиеся части массивов ещ егзгз будут сдвинуты вверх на одну группу запаздывания. Величины в последних группах будут заменены нулями.
Счет переноса энергии на одном шаге по времени
Вычисления, связанные с моделированием переноса, на одном шаге по времени минимальны:
вычисление энергии, излучаемой за шаг по времени сеточными источниками, — на основе значений плотности Jj" и uT)STl с предыдущего шага, (79);
распределение (умножение на коэффициенты переноса) проникающей части энергии, размещение вычисленных импульсов энергии в массиве ещ распределение выведенной веществом части энергии, размещение вычисленных импульсов энергии в массиве ew; перенос первой группы запаздывания массива еи в массив Jy, сдвиг остальной информации на одну группу запаздывания;
счет граничных условий, вычисление нового значения плотности J^" излучаемой энергии для граничных источников, (9);
перенос первой группы запаздывания массива erw в массив uaiSUl сдвиг остальной информации на одну группу запаздывания;
счет поглощенной иаТ1 и рассеянной энергии и8Нч (84);
решение итерациями Ньютона уравнения баланса энергии для вещества, счет новых значений коэффициентов взаимодействия излучения и вещества, вычисление нового значения плотности ит 8П для внутренних источников, (35)-(37).’
Обсуждение результатов
Нам осталось определить схему вычисления базовых коэффициентов, Pip Pip Pin, Pin, Ртр Ртр Ртпч Ртп* На ИХ ОСНОВе рассчитываются (п. 4.6) коэффициенты переноса, рг-к, пру, ррк; Pink’ ^гпкч Pink’ Pmjk1 ^mjki Ртркл Pmnk’ ^тпкч Ртпк'
Обе группы коэффициентов, базовых и переноса, позволяют описать (п. 4.3) распределение импульса сеточных источников: проникающей части энергии — по граничной поверхности, выведенной веществом части энергии — внутри объема вещества. Описание переноса с помощью коэффициентов согласовано с шагами по времени.
Кроме описания переноса мы получаем равенства (85), (86) для сумм коэффициентов распределения.
Сумма точных значений должна быть равна единице. Сумма приближенных значений отличается от единицы.
Величина суммы приближенных значений, с одной стороны, служит критерием точности счета, с другой стороны, — основой для нормирования коэффициентов распределения. Деленные на эту величину коэффициенты дают возможность выполнения счета переноса по явной схеме с полным балансом энергии для каждой сеточной ячейки в отдельности.
Это совершенно новый уровень балансности схемы счета переноса и контроля точности счета, по сравнению со счетом с помощью итераций «по правой части» на основе интегро-дифференциального уравнения или интегрального уравнения Пайерлса.
Коэффициенты ослабления, базовые и переноса, связывают, в представленных выше формулах, все счетные ячейки-источники со всеми счетными ячейками-целями. В результате схема свободна от ошибок, связанных с лучевыми эффектами счета.
Схема счета коэффициентов ослабления, к которой мы перейдем ниже, является численной аппроксимацией физической корпускулярной модели переноса. В то же время, схема вычисления базовых коэффициентов — схема приближенного вычисления определенных интегралов по поверхности или по телесному углу.
Эти интегралы определены, в простейшем случае, законом Ламберта описания распределения для граничных источников, и законом Бугера-Ламберта.
Глава 5
ВЫЧИСЛЕНИЕ КОЭФФИЦИЕНТОВОСЛАБЛЕНИЯ
Коэффициенты — двукратные определенныеинтегралы
В ограниченной области пространства V с граничной поверхностью S численно решается нестационарная краевая задача переноса излучения с помощью сети ячеек: внутренних ячеек Vm, V = J^m, и граничных ячеек Si, S = (J^. Теми же буквами
обозначены объем Vm внутренних ячеек и площадь Si граничных ячеек. Счет ведется шагами по времени, шаг постоянен. Все заданные и рассчитываемые функции кусочно-постоянны по пространственным переменным и времени: они постоянны в пределах одной ячейки и в пределах одного шага по времени. Но могут иметь другие значения в пределах другой ячейки и в пределах другого шага по времени.
Счетные центры ячеек Pi, Рт — место сеточных точечных граничных и внутренних источников лучистой энергии, описывающих излучательную способность ячеек, Pi е Si, Рт Е Vm.
В такой сеточной постановке краевая задача переноса лучистой энергии является квазианалитической [8], она сводится к вычислению определенных интегралов по телесному углу и к решению итерациями Ньютона нелинейных граничных уравнений и уравнений баланса энергии вещества. Поэтому краевая задача в сеточной постановке может быть решена с гарантированной контролируемой точностью.
Схема вычисления базовых коэффициентов Д^, р^, р^, Pink, Pink', Pmjk, Pmjk, Pmjk, Pmnk, Pmnk ЯВЛЯСТСЯ ЧИСЛвННОЙ Трехмерной аппроксимацией физической корпускулярной модели переноса импульсов, равных единице. Вместе с точкой-источником рассматривается пучок расходящихся лучей, ведущих в центры элементов разбиения граничных ячеек. Дополнительное разбиение ячеек на элементы — средство численного интегрирования
по поверхности ячеек — переход от интеграла к допредельной частичной сумме.
Схема вычисления базовых коэффициентов — схема прибли
женного вычисления определенных интегралов по поверхности граничных ячеек или по боковой поверхности внутренних ячеек,
частичаститочнее, — интегралов по телесному углу.
части
части
Для коэффициентов распределения проникающей
энергии:Pijk =энергии:Pijk =Pijk —Pijk —Sjk(Pi)j exp(-TP1Q)KPiQdSQ,Sjk(Pi)j exp(-TP1Q)KPiQdSQ,^^dSQ,ppp.o^^dSQ,ppp.oSjk(Pt)Pijk = -P [ PPiQexp^-Tp^Kp^dSQ, Pi IkJSik^Ppj KWSq. K,p„.Qd-.Q=3P-Sjk(Pm)^PPmQPmjk= j exp(-TPmQ)KpmQdSQ,Sjk(Pm)Pmjk=P— [ PPmQeXP(-TPmQWpmQdSQ, Pmj кJSjk(Bm)Sjk(Pt)Pijk = -P [ PPiQexp^-Tp^Kp^dSQ, Pi IkJSik^Ppj KWSq. K,p„.Qd-.Q=3P-Sjk(Pm)^PPmQPmjk= j exp(-TPmQ)KpmQdSQ,Sjk(Pm)Pmjk=P— [ PPmQeXP(-TPmQWpmQdSQ, Pmj кJSjk(Bm)dSQ,dSQ,(87)(87)
энергии:
Pijk =
энергии:
Pijk =
Pijk —
Pijk —
Sjk(Pi)
j exp(-TP1Q)KPiQdSQ,
Sjk(Pi)
j exp(-TP1Q)KPiQdSQ,
^^dSQ,ppp.o
^^dSQ,ppp.o
Sjk(Pt)
Pijk = -P [ PPiQexp^-Tp^Kp^dSQ, Pi IkJ
Sik^Pp
j KWSq. K,p„.Qd-.Q=3P-
Sjk(Pm)^PPmQ
Pmjk= j exp(-TPmQ)KpmQdSQ,
Sjk(Pm)
Pmjk=P— [ PPmQeXP(-TPmQWpmQdSQ, Pmj кJ
Sjk(Bm)
Sjk(Pt)
Pijk = -P [ PPiQexp^-Tp^Kp^dSQ, Pi IkJ
Sik^Pp
j KWSq. K,p„.Qd-.Q=3P-
Sjk(Pm)^PPmQ
Pmjk= j exp(-TPmQ)KpmQdSQ,
Sjk(Pm)
Pmjk=P— [ PPmQeXP(-TPmQWpmQdSQ, Pmj кJ
Sjk(Bm)
dSQ,
dSQ,
(87)
(87)
Для коэффициентов распределения выведенной веществом части энергии излучения:
Гпк = j exp^-Tp^Kp^dSQ,
Гпк= j exp^-Tp^Kp^dSQ,
ptnk= j PP^exp^-Tp^Kp^dSg,
P?nk= j Pp^exp^-Tp^Kp.QdSq,
Snk(Pi)
Pink ^ink ^ink’
Pmk = O,b^nkl0tnk^ PlnklSinkY uink = 0M3tnk"V ЗГпкУ (88) 3mnk=|exp^-Tp^KPmQdSq,
3mnk=|^-TPmQXPmQdSQ,
Pmnk =|№mQ ^M-TPmQWPmQ ^Q,
S^(Pm)
Pmnk =|№mQ ^Vk-TPmQWPmQ ^Q^
Pmnk Pmnk Pmnk1
Pmnk ^^^Pmnk/Pmnk + Pmnk! PmnkP ^mnk ^P^Pmnk + PmnkP Коэффициенты pijk, pijk, pijk, piukl Pink описывают распределение импульса излучения для граничного сеточного источника:
— Pijk — доля углового распределения энергии излучения граничного сеточного источника Р^,
угловой коэффициент переноса проникающей части энергии в сторону граничной ячейки-цели Sjk,
— Pijk — доля полного распределения энергии излучения граничного сеточного источника
коэффициент переноса проникающей части энергии в сторону граничной ячейки-цели Sjk с учетом ослабления излучения веществом,
— Pijk — среднее расстояние от граничного сеточного источника Pt до граничной ячейки Sjk,
— Pink — суммарная доля углового распределения энергии для лучей граничного сеточного источника Р^, пересекающих внутреннюю ячейку Vnkl
угловой коэффициент переноса выведенной веществом части энергии в сторону ячейки-цели Vukl
— Pink — доля полного распределения энергии излучения граничного сеточного источника Р^
коэффициент переноса (выведенной части энергии) в сторону внутренней ячейки Vuk с учетом ослабления излучения веществом этой ячейки,
— Pink — среднее расстояние от граничного сеточного источника Pt до внутренней ячейки Vnk.
Коэффициенты pmjk, pmjk, pmnk, pmnk описывают распределение импульса излучения для внутреннего сеточного источника Рт.
Из (71), (72), (87), (88) следуют контрольные соотношения (85), (86) для коэффициентов распределения энергии излучения граничного, p^kl Рть и внутреннего, pmjkl ртпк, сеточных источников.
Для угловых коэффициентов Д^, 3mjk контрольные соотношения имеют вид:
У2^к = | KpiQdSQ=l1 ^23mjk= | KPmQdSQ=L jk S(pp 3k s^
(89) что проверяется непосредственным интегрированием, см. введение.
Функция Кр^ описывает угловое распределение для граничного точечного источника как произведение двух вероятностей: вероятности излучиться в фиксированном направлении PtQ для поверхности, излучающей по закону Ламберта (пропорционально рЗь и вероятности излучиться внутрь телесного угла (1ш9, duQ = pQdsQfp^Q.
После перехода в местную цилиндрическую систему координат, из (89) непосредственным интегрированием получаем величины нормирующих множителей, Кр^ = costM^q/tt, dujQ = smddddp, ^'pmQ = ^q/(4tt).
Функция KpmQ описывает угловое распределение для изотропного внутреннего точечного источника как вероятность излучиться внутрь телесного угла.
Фиксируя номер сеточного источника, г или т, тем самым фиксируем счетный центр — точечный источник Рг или Рт.
Можно было бы отдельно вычислять коэффициенты распределения проникающей и выведенной веществом частей энергии для каждой граничной и внутренней ячейки, следуя формулам (87), (88).
Например, можно вычислить по каждому дискретному лучу проникающие части импульса на входе и выходе для освещенной и теневой сторон поверхности пересекаемой ячейки отдельно, и потом найти разность — выведенную веществом энергию. При этом необходимо для определения оптической толщины пути до входа, повторять вычисления координат точек пересечения лучей с границами других предшествовавших внутренних ячеек.
Введем понятие согласованного счета базовых коэффициентов:
для граничных сеточных источников — коэффициентов Pijki Pijk, Pijk описания распределения проникающей части энергии по граничным ячейкам и коэффициентов Pinp, Piny описания распределения выведенной веществом энергии;
для внутренних сеточных источников — коэффициентов 3mjb Pmjki Pmjk описания распределения проникающей части энергии по граничным ячейкам и коэффициентов Ртпк, Pmnk описания распределения выведенной веществом энергии.
Для численного интегрирования рассмотрим некоторое дополнительное разбиение граничных ячеек-целей Sjp на элементы Sjpq, Sjp = IJS^g, с выделением в каждом из них счетного q
центра Qq — для определения счетного луча PQq и вычисления подынтегральной функции.
При выборе дискретных средств для приближенного вычисления определенных интегралов должно быть выполнено требование стремления к нулю диаметра областей разбиения. Оно будет выполнено, если разбиение граничной поверхности и внутренних дискретных объемов проводится — даже для двумерных расчетов — по всем трем координатам. Например, для двумерных расчетов — по цилиндрическим координатам, включая угол вращения.
В двумерном случае с учетом осевой симметрии вращения при дополнительном разбиении ячеек-целей по углу вращения вычисления производятся в диапазоне изменения угла вращения 0 < (/> < тг.
Рассмотрим базовое множество лучей (гл. 3) — множество лучей PQq, ведущих в счетные центры Qq элементов разбиения граничной поверхности. Базовое множество лучей определяет базовое разбиение dojq телесного угла для точки-источника Р.
Здесь dwq = dsQ^Qjp2PQq.
Согласованным счетом всех интегралов для половины телесного угла (2л) — для граничных сеточных источников или для полного телесного угла (4л) — для внутренних сеточных источников будем называть их вычисление с использованием только базового множества лучей и базового разбиения телесного угла. При этом используется «наведенное» разбиение поверхностей S\ S~ на элементы dsQ^, dsQ-, dwq = duQ^ = duQ-, см. рис. 1.
По каждому лучу PQq базового множества вычисления проводятся один раз: от точки-источника Р до точки Qq на границе области.
Последовательно вычисляются координаты точек входа и выхода для пересекаемых лучом внутренних ячеек, оптическая толщина пути и ослабление.
На основе закона Бугера-Ламберта определяется энергия, выведенная веществом этих ячеек. Она становится слагаемым при вычислении коэффициентов распределения выведенной энергии для пересекаемых ячеек.
Наконец, оставшаяся энергия (от введенной единицы) становится слагаемым при вычислении коэффициентов распределения величины Jp(Q,t + Ppq/с) граничной энергии — проникающее излучение.
Проводим лучи PiQq. Когда точка Qq пробегает центры разбиения граничной поверхности, доли телесного угла, опирающиеся на элементы разбиения, составляют половину телесного угла.
В двумерном случае поверхности Sjk — части поверхности вращения. Разбиение их на элементы Sj^q определяются двумя целыми числами: числом ка разбиения с равным шагом по длине образующей и числом к- равномерного разбиения по углу вращения, так же с равным шагом. Каждая граничная ячейка-цель разбивается на ка х к^ элементов численного интегрирования.
Приближенные значения коэффициентов удовлетворяют контрольным соотношениям с некоторой ошибкой. Отличие сумм коэффициентов распределения от единицы служит характеристикой точности счета. Кроме того, вычисленные суммы коэффициентов служат основой для поправочной нормировки приближенных значений коэффициентов. Нормировка значений коэффициентов не повышает уровень точности счета, но приводит к точному балансу энергии, недоступному для предшествовавших численных методов. Балансность нового метода локальна для каждого сеточного источника. В предшествовавших методиках с помощью итераций совместного решения полной системы приближенных уравнений достигается приемлемая точность выполнения суммарного баланса.
Анализ проверки выполнения контрольных соотношений (до нормировки коэффициентов) для тестовых задач дал возможность сделать некоторые выводы о свойствах новой численной схемы.
Основные неточности счета связаны с описанием углового распределения.
Контрольные соотношения (85), (86) для коэффициентов распределения энергии излучения граничного, р^у, Рть и внутреннего, pmjk, Pmnk, сеточных источников выполняются с той же точностью, что и контрольные соотношения (89) для коэффициентов Р^у, Pmjk углового распределения. Вычисления ослабления энергии при движении вдоль луча выполняются с точностью производимых геометрических вычислений и не изменяют число точных значащих цифр результатов счета.
Первым выполняется счет угловых коэффициентов Д^, 3mjk и проверяется выполнение контрольного соотношения (89). Можно установить режим пересчета углового распределения до достижения заданной точности выполнения контрольных соотношений. Эти вычисления выполняются итерациями с ростом чисел kai k^, определяющих разбиение ячеек граничной поверхности для плохо сосчитанных источников.
Сами угловые коэффициенты З^у, 3mjk ПРИ счете задачи переноса энергии не нужны, используются коэффициенты, нормированные с учетом контрольных соотношений.
С помощью установленного разбиения ka х fc^ выполняется численное интегрирование для расчета всех коэффициентов распределения энергии излучения для фиксированного сеточного источника.
Алгоритм пересечения лучом внутренних ячеек
Закрепим сеточные числа:
г или т, определяющие граничную или внутреннюю ячейку, сеточный источник излучения Р^ или Рт,
у, определяющее граничную ячейку-цель Sj,
к, определяющее долю Sjk разбиения ячейки-цели по углу вращения — для вычисления среднего расстояния и учета запаздывания,
q, определяющее площадь Sjkq и ее счетный центр Qq при разбиении доли Sjk на ка х к^ элементов численного интегрирования.
Тем самым, мы фиксировали базовый счетный луч PtQq ИЛИ PmQq-
Координаты точек Р^ Рт — сеточных источников излучения — рассчитываются на основе декартовых координат узлов сети Х(пп), ¥(тш) на стадии подготовки расчета, пи — номер узла сети.
Кроме массивов Х(тш), Ч(пи), содержащих всю геометрическую сеточную информацию о задаче, на стадии подготовки
расчета формируется информация о топологических функциях узлов и образующих отрезков, а также описание контуров внутренних ячеек. Все вместе эти массивы позволяют обеспечить расчет коэффициентов в условиях использования широкого класса двумерных сеток.
Массивы NUN(no), NUR(no) сообщают для каждого номера по образующего отрезка в полуплоскости р = 0, р — угол вращения, номера узлов — точек начала и конца отрезка. Первыми описываются граничные образующие отрезки, 1 < по < кд. Для описания вектора внутренней нормали для граничных отрезков важна векторная форма описания: предполагается, что расчетная область расположена справа от образующих граничных отрезков.
Массивы Ml (по), М2(по) сообщают номера внутренних ячеек, разделяемых каждым образующим отрезком с номером по. Для граничных отрезков один из номеров примыкающих ячеек равен нулю.
Контуры внутренних ячеек описываются двумя массивами NOM, АОМ(ш).
Массив NOM неоднороден, для каждой ячейки с номером т сообщаются:
число ко образующих отрезков описания контура ячейки в полуплоскости р = 0, например, ко = 4 — для четырехгранных ячеек,
затем перечисляется ко чисел — номера образующих отрезков.
Массив АОМ позволяет считывать в массиве NOM информацию о контуре внутренней ячейки с номером т, она начинается с адреса АОМ(т).
Остальные массивы связаны с описанием постановки краевой задачи: описанием параметров вещества и граничных условий, заданием начальных значений рассчитываемых величин.
Все массивы одномерны.
С помощью массивов координат точек Р^ Рт и массивов NUN, NUR на основе теоремы косинусов для треугольников вычисляется коэффициент углового распределения для луча, Кр^ ИЛИ KPmQ'
Далее рассчитывается перенос порции энергии, Кр^ или KpmQ, вдоль луча.
Первой пересекается внутренняя ячейка с номером т, примыкающая к граничному отрезку с номером г, или внутренняя ячейка сеточного источника Рт.
Точка Pi или Рт является точкой Ро входа луча PiQq или PmQq в первую пересекаемую внутреннюю ячейку. Вычисляются координаты точки выхода Рх луча из этой ячейки.
Для этого вычисляются точки пересечения луча с каждой из граней — поверхностей вращения образующих отрезков, составляющих контур внутренней с номером т. Используются массивы AOM(m), NOM, NUN, NUR, X, У.
Точкой выхода является точка пересечения, ближайшая к точке входа.
Вычисляется длина пути PiPo и Рг?х (или РтРо и РтРх) от сеточного источника до точек входа и выхода, оптическая длина пути PqPx по пересекаемой ячейке и ослабление ехр(—тр.р0), ехр(—rp.pj на пути до точек входа и точек выхода. Вычисляется адрес в массивах базовых коэффициентов на основе чисел номеров источника и пересекаемых долей разбиения внутренних ячеек по углу вращения.
Все вычисленные величины не переносятся сразу в массивы вычисляемых коэффициентов, а запоминаются сначала в рабочих массивах, рассчитанных на число внутренних ячеек и их долей разбиения, пересекаемых одним лучом.
После занесения вычисленных величин для первой пересекаемой ячейки переходим к следующей. Номер следующей пересекаемой ячейки определяется с помощью массивов Ml (по), М2(по) по номеру по отрезка образующей для грани, пересекаемой в точке выхода.
Точка выхода становится точкой входа в соседнюю пересекаемую внутреннюю ячейку, точнее в соседнюю долю разбиения пересекаемой ячейки. И так до исчерпания, пока точка выхода не совпадет с граничной точкой конца луча.
Лучевой алгоритм видимости
Пока мы считали граничную поверхность выпуклой.
В случае невыпуклости граничной поверхности встает вопрос о вычислении области видимости S^P) на границе для точечного источника Р.
Некоторые лучи могут быть направлены в точки Qq границы, затененные самой граничной поверхностью. Двигаясь вдоль такого луча от сеточного источника, мы придем в точку Рх выхода на границу, Рх 6 S, не совпадающую с точкой конца луча, Pi Qq
В этом случае луч выходит из расчетной области, он не актуален. Вычисленные результаты не переносятся из рабочих массивов в массивы рассчитываемых базовых коэффициентов.
Край области видимости при этом определяется с первым порядком точности О(А^), где А92 — шаг по углу вращения для элемента разбиения S^q со счетным центром Qq при численном интегрировании. Вопрос о включении слагаемых — промежуточных величин — в результат, определяется для центрального луча PiQq или PmQq, а относится к телесному углу, опирающемуся на площадь элемента разбиения As = Аа х А99, Аа — шаг вдоль образующего отрезка.
При решении двумерных задач можно повысить порядок точности описания области видимости 8(Р\
Очередной шаг Аа вдоль образующей определяет «поясок» — соответствующую часть поверхности вращения Sjk- Геометрическим местом точек Qqi относящихся к «пояску», является окружность — центральная линия «пояска». На этой окружности с помощью алгоритмов, подобных описанным в [15, 21], можно аналитически вычислить границы дуг видимости для точки- источника, Рг или Рт. Разбиение «пояска» по углу вращения можно выполнить с учетом этих границ.
Глава 6
ЧИСЛЕННОЕ РЕШЕНИЕ ТРЕХМЕРНОЙ ЗАДАЧИЭНЕРГОПЕРЕНОСА ИЗЛУЧЕНИЕМ ВНУТРИПОЛОСТИ СЛОЖНОЙ ФОРМЫ С УЧЕТОМ
ЗАПАЗДЫВАНИЯ И ЗАТЕНЕНИЯ
Аппроксимация граничных поверхностей.Граничные условия
Остановимся на решении двумерных и трехмерных задач нестационарного переноса лучистой энергии в вакуумном приближении. Будем считать, что внутри объема расчетной области нет заданных источников излучения, а взаимодействием излучения и вещества среды можно пренебречь, xaiS = 0.
В этих условиях отсутствуют и наведенные источники излучения внутри объема, uT)S = 0. Есть только граничные источники излучения, обменивающиеся энергией через полость сложной формы — с учетом запаздывания и затенения. Геометрия расчетной области не меняется в ходе расчета.
В двумерном осесимметричном случае граница аппроксимируется участками поверхностей вращения. Задача вычисления точек пересечения границы и лучей переноса энергии сводится к решению квадратных уравнений.
В [15, 18] для аппроксимации границ в полуплоскости у? = 0 используются ломаные линии из отрезков прямых и дуг окружностей с центром на оси вращения. Отрезки и дуги служат образующими участков поверхностей вращения: конусов, цилиндров, сфер и «дисков» — колец в плоскостях, перпендикулярных оси вращения. Для более аккуратного учета эффекта запаздывания используется дополнительное разбиение поверхностей вращения по углу вращения.
Для определения образующего отрезка прямой или дуги сферы в полуплоскости у? = 0 достаточно указания координат двух точек: точки начала и точки конца. Отрезок или дуга становятся при этом направленными, векторными объектами. Для определения расчетной области при описании аппроксимирующих гра-
104 Гл. 6. Численное решение трехмерной задачи энергопереноса ницу ломаных линий должно соблюдаться выбранное правило обхода области, левое или правое.
При решении трехмерных задач переноса лучистой энергии в [18] граница описывается с помощью основной и вспомогательных систем декартовых координат. К описанным выше участкам поверхностей вращения с заданными пределами изменения угла вращения добавлены участки координатных плоскостей. Указывается число разбиений этих участков на счетные ячейки. Ячейки нумеруются в порядке описания аналитических участков граничных поверхностей.
Соблюдаются правила обхода расчетной области. При описании стыковки участков, определенных в разных системах координат, возникают трудности. Требуется избежать использования кривых четвертого порядка.
При решении двумерных осесимметричных задач движущиеся границы можно аппроксимировать с помощью поверхностей вращения отрезков прямых. Координаты узлов могут, например, вычисляться в ходе совместного расчета движения. Особенно удобны расчеты с использованием лагранжевых сеток счетных ячеек. Координаты узлов образующих отрезков являются непрерывными функциями времени, кроме моментов времени, обусловленных необходимыми перестроениями сети ячеек.
При решении трехмерных задач поверхности вращения, использованные для аппроксимации начальной геометрии расчетной области, могут не сохраняться в процессе движения. Встает вопрос об аппроксимации поверхности граней. Например, плоскостями треугольников, построенных с помощью точек-узлов грани и точки их центра тяжести [36]. Если в программе расчета движения вещества среды, кроме вычисления координат узлов сети, предусмотреть определение координат центра грани с учетом стрелы прогиба этой грани, можно повысить порядок точности схемы счета. Описание трехмерных поверхностей с помощью треугольных элементов используется в комплексах построения трехмерной графики.
В [18] речь идет о решении задач с неподвижной геометрией расчетной области.
Граничные условия описывают (см. [4, 7]) заданные источники энергии, либо баланс энергии: пришедшая на поверхность граничной ячейки энергия частично отражается обратно внутрь расчетной области, частично отводится в окружающее пространство,
J~(P,t) = J+(F,£) + g(F,£), Ре Si.
Предполагается также заданной, аналитически или таблично, функция f(ptp) углового распределения излучения граничного источника; pip — косинус угла наклона направления переноса к вектору внутренней нормали в точке Р. Для изотропного излучения f(ptp) = Hpfiv (распределение Ламберта).
Из интегро-дифференциального уравнения в случае вакуумного уравнения следует, что при переносе вдоль луча интенсивность излучения сохраняется.
Потоковые модели при решении задач в вакуумном приближении могут обеспечить сохранение интенсивности лишь на некоторых выделенных направлениях — с помощью закрепления этих направлений в качестве сеточных и с помощью удачного подбора сеток, использующих особенности геометрии расчетной области. Успех расчета в этом случае является следствием искусства исполнителя расчета. Связать все граничные ячейки лучами-характеристиками в рамках потоковой модели практически невозможно из-за использования единого раскроя телесного угла в счетных центрах всех ячеек.
Для решения задач в вакуумном приближении удобны интегральные уравнения переноса излучения.
Аппроксимация интегрального уравненияПайерлса
Особенности гладкости решения уравнения переноса в вакуумном приближении есть особенности гладкости одного из двух слагаемых интегрального уравнения переноса в общем случае, [13]. Второе слагаемое — энергия излучения внутренних источников — отсутствует.
Для граничной точки-цели Q интегральное уравнение Пайерлса (61) в этом случае имеет вид:
JpQ.t^ | P(P,t-ppQ/c)KpQdSp, (90)
S(Q)
Kpq = ^yP pp = cos^p,PQ\ pq = cos(iiq,QP).
На единицу площади в окрестности граничной точки-цели Q в момент времени t переносится энергия излучения с видимой части граничной поверхности с учетом запаздывания. В случае невыпуклости граничной поверхности актуальны лишь лучи PQ, ведущие в точку-цель Q без пересечения границы во внутренних точках пути PQ. По разным направлениям — разные сдвиги
по времени ppq/c. Необходимо помнить «историю» — функции состояния источников в предыдущие моменты времени.
Характерным становится диапазон Dt времени, соответствующий времени прохождения светом диаметра расчетной области. Используется предположение о квазистационарности границ расчетной области — о неподвижности границ на время Dt — в связи с большим различием по величине между скоростью света и скоростью движения вещества.
В этом предположении уравнение (90) в вакуумном приближении аппроксимируется алгебраическим уравнением переноса лучистой энергии (91), см. ниже.
С учетом разбиения граничной поверхности на ячейки и сеточным кусочно-постоянным описанием физических величин, уравнение (90) имеет вид:
J (Q,t) = ^2 | J^(P,t - ppQ/c)KpQdSp, г S^Q)J (Q,t) = ^2 | J^(P,t - ppQ/c)KpQdSp, г S^Q)ASqSdQ)ASqSdQ)QeSpPeSvQeSpPeSvKpq dSp,Kpq dSp,
J (Q,t) = ^2 | J^(P,t - ppQ/c)KpQdSp, г S^Q)
J (Q,t) = ^2 | J^(P,t - ppQ/c)KpQdSp, г S^Q)
ASq
SdQ)
ASq
SdQ)
QeSpPeSv
QeSpPeSv
Kpq dSp,
Kpq dSp,
Здесь t — шаг по времени, n — номер текущего шага по времени, tn — момент времени, соответствующий текущему шагу. Выполнено дополнительное интегрирование уравнения по поверхности ячейки-цели Sp интегрирование проводится по той части граничной поверхности, для которой видят друг друга сеточные центры разбиений ячейки-цели и ячейки-источника. Энергия, приходящая в граничную ячейку-цель за время г, складывается из долей энергии, излучаемой источниками в сторону этой ячейки,
г
сьц = |j KpqdSpdSq,
S3(P)S4(Q)
(91)(91)Pij = — [[ PPQKpQdSpdSQ,
(91)
(91)
аг3
p W = 52 Pij Г Г - pijfdSiT,
=
Область видимости описывается симметричным образом: интегрирование ведется по участкам ячеек Si, Sj, Q Е Sj, Р Е Si, находящимся во взаимной видимости. В силу симметрии подынтегральной функции aij = aji. Коэффициенты aij называются коэффициентами взаимной видимости двух граничных ячеек.
Из (91) при изотермии — для системы, прогретой до стационарной температуры, — следует равенство ^aij = Sj, ^ciij =
г 3
= Si. Эти же равенства могут быть установлены непосредственным интегрированием. Эти равенства служат контрольными соотношениями, точность их выполнения характеризует точность расчета.
Доли pij называются коэффициентами углового распределения энергии ячейки-источника Si в сторону ячеек-целей Sj, или угловыми коэффициентами, ^Pij = 1.
з
Если шаг по времени достаточно мал и ст < р^, то справа в (91) все величины вычислены на предыдущих шагах, следовательно, вычисления носят явный характер. Иначе часть величин или все величины определяются через неизвестные величины и вычисления неявно. Для общности предполагается последний вариант. Совместно решается полная система уравнений, для ее решения используются итерации. Начальными значениями для всех неизвестных величин служат их значения с предыдущего шага. Шаг т меняют в ходе расчета с целью повышения экономичности расчета, но с внесением ошибок первого порядка точности О(т) при каждой смене шага.
С помощью нормирования угловых коэффициентов, р'- = = р^/^2Pij, можно достичь точного выполнения соотношения
^2 Pij = 1- При использовании явного варианта счета с корректи- г
рованными коэффициентами р'- вместо pij, при постоянстве шага по времени и при постоянстве по времени геометрии границ, счет будет балансным.
В общем случае, при совместном решении полной системы аппроксимирующей уравнений с помощью глобальных итераций, счет не обладает балансностью. Источники востребуются в разные моменты времени с разными угловыми коэффициентами, вычисленными для разных состояний геометрии расчетной области.
Аппроксимация интегрального уравненияв новой форме
Новая форма интегрального уравнения (71) для граничного точечного источника Р в случае вакуумного приближения имеет вид:
Р(Р,1р8рт = | JpP.t + ppq/cpSQT. (92)
S(F)
Импульс энергии распределился по граничной поверхности и по времени, с учетом запаздывания. Здесь Jp(Q,£ + Ppq(c^ — доля точечного источника Р в плотности одностороннего потока энергии на границе области;
Jp(Q, t + ppQ/c)dSQ = J^(P, t^dSpKpQdSQ,KPQdSQ = t^dsq = ^Л±ЛЛ
— функция углового распределения энергии точечного источника Р.
Сделаем предположение, обычное для сеточных методов численного описания переноса излучения: импульс энергии излучения ячейки St высвечивается из точечного сеточного источника Pi, счетного центра этой ячейки.
Подготовленный предыдущим изложением читатель должен рассматривать предыдущие формулы, начиная с (92), как формулы прямого действия. С появлением ЭВМ появился круг квазианалитических средств, использование которых выполнимо с гарантированной заказанной точностью, [9]. Следующая часть изложения является повторением сказанного для сеточного дискретного описания функций с использованием численного интегрирования.
С учетом разбиения на ячейки граничной поверхности и сеточным кусочно-постоянным описанием физических величин уравнение (92) примет вид:
pPPdt = 'p | J^qP + pp^/cpsqdt =
= E J^HSidt | KPiQ dsQ dt, (93)
pPpdt = 'ppijJppSidt, ^Pij = I- jj
Здесь слева — величина импульса энергии, излучаемой граничным источником Pt за шаг по времени; Jpg — вклад этого источника в энергию, приходящую на элемент dsg граничной ячейки-цели Sj; р^ — угловой коэффициент, доля углового распределения импульса энергии,
f Pp1Pq л Г 7
Pij = J pp^dsQ=
d^PiQ — телесный угол, под которым из точки-источника Pi виден элемент dsg граничной ячейки-цели Sj.
При аппроксимации уравнения Пайерлса в вакуумном приближении угловые — четырехкратные интегралы, по поверхности Si ячейки-источника и по поверхности Sj ячейки-цели. В нашем случае — при аппроксимации интегрального уравнения в новой форме — коэффициенты рр распределения — двукратные интегралы.
Интегрирование выполняется по области видимости поверхности ячейки-цели Sj^Pi) для точечного источника излучения Pi, счетного центра граничной ячейки Si.
На современных ЭВМ вычисление определенного двукратного интеграла может быть выполнено с квазианалитической [9] точностью — с высокой контролируемой точностью.
Схема счета (92), (93) переноса в вакуумном приближении на основе интегрального уравнения (71)-(75) в новой форме есть часть схемы счета переноса лучистой энергии на той же основе в общем случае, с учетом взаимодействия излучения с вещества среды.
За шаг по времени т сеточный граничный источник Si излучает энергию Е{, ЕТ = J^SiT.
Для определения доли Е^ этой энергии, переносимой внутри полости с учетом запаздывания к цели Sj, воспользуемся фрагментом нового интегрального уравнения, опуская описание ослабления и вынося за знак интеграла среднее значение излучаемой энергии:
Eij(tpij/с) = Ei рр, pij = | Kp^ dSg, Q 6 Sj.
Чтобы «доставить» до цели вычисленные доли энергии Е^ используется массив ер [18], в котором переносимая энергия распределена по группам запаздывания, а внутри группы — по целям.
Особенности решения трехмерных задач
При вычислении трехмерных коэффициентов видимости отметим некоторые моменты, связанные с численным интегрированием по области видимости Sj(Pi) на поверхности ячейки-цели для фиксированного источника Pi с учетом затенения невыпуклой граничной поверхностью сложной формы.
Определены координаты точек начала и конца для луча источник-цель PiQjq. Точка Pi — излучающий центр счетной ячейки-источника Si. Точка Qjq — центр очередного элемента Sjq разбиения поверхности счетной ячейки-цели Sj, Sj = U*SjQ, — q
для выполнения численного интегрирования (75). Нужно проверить, что луч PiQjq не имеет внутренних пересечений с границей — в случае невыпуклости границ.
При счете с неподвижной геометрией границ можно для сокращения вычислений особо выделить участки Sm аналитического описания границы, мешающие видимости. В один «мешающий» участок Sm можно объединить несколько граничных ячеек, Sm = (jSJm, с единым аналитическим описанием (части
Зт
одного конуса, например).
Чтобы не выходить за пределы решения квадратного уравнения, приходится переводить координаты точек Pi, Qjq в систему координат очередного участка Sm, мешающего видимости. Участок Sm — поверхность вращения второго порядка, аналитически определенная в некоторой декартовой системе координат, основной или вспомогательной. В любой из этих систем координат луч PiQjq определяется уравнением прямой, проходящей через точки Pi, Qjq. Координаты точек Pi, Qjq переводятся в систему координат описания участка Sm, мешающего видимости.
Точка Pi «видит» точку Qjq, если отрезок PiQjq не имеет внутренних пересечений с границей. Будем считать, что в этом случае точка Pi «видит» весь элемент разбиения Sjq, Sjq С Sj^Pi). Область видимости Sj^Pi) при этом с первым порядком точности описывается совокупностью видимых элементов, Sj^Pi) = |J Sjq(Pi).
В [15, 21] при решении двумерных задач используется другой алгоритм описания области видимости. Там использован фактор соосности поверхностей вращения для ячейки-цели и для «мешающей» поверхности.
Изложим аналогичный алгоритм в применении к решению трехмерных задач. Пусть граничные поверхности аппроксимированы треугольными элементами. S = |J Аа, S3 = |J Aja.
Можно с помощью проекций «мешающих» треугольных элементов на плоскость фиксированного треугольного элемента определить точную область видимости. Это будет один или несколько многоугольников сложной формы, вообще говоря, невыпуклых. Нужно будет преодолеть трудности их описания и разбиения на элементы для численного интегрирования.
В соответствии с [15] разобьем поверхность фиксированного треугольного элемента Aja на полосы — треугольники с общей вершиной. Проведем внутри каждого из них средние линии — диагонали углов при общей вершине. С помощью проекций «мешающих» треугольных элементов на плоскость фиксированного треугольного элемента определим точные границы отрезков видимости на всех средних линиях. После этого разбиения на элементы численного интегрирования полос будем производить с учетом вычисленных точных границ отрезков видимости для средних линий. При этом площадь совокупности видимых элементов Sjq(Pi) отличается от точной области видимости на величину второго порядка малости 0(5^), Sj(Pi) = |J Sjq(Pi) ± 0^\ где 6q — наибольший диаметр разбиения поверхности счетной ячейки S3 для выполнения численного интегрирования (91). Площадь видимой части ячейки-цели вычислена по методу средних прямоугольников.
Очень важны предварительные проверки, исключающие рассмотрение тех или иных треугольников в качестве «мешающих». Например, расположенных по какой либо координате слева или справа от минимума-максимума значений этой координаты для луча PiQjq. В качестве «мешающих» видимости можно рассматривать лишь треугольники с тупым углом наклона вектора нормали к вектору PiQjq.
Необходимость в использовании алгоритма вычисления области видимости с повышенной точностью особенно ощутима в расчетах сложных затенений. Примером могут служить расчеты полостей, перегороженных шторками. Источник излучения и облучаемая поверхность находятся по разные стороны от шторки [15]. В этом случае эффективно задание поверхностей-«фильтров», пересечение которых необходимо для лучей, реализующих перенос лучистой энергии между подобластями полости. В трехмерном расчете, описанном в [15], две осесимметричные системы были разделены непрозрачной пере-
городкой с отверстием квадратного сечения. Сечение отверстия было описано как поверхность-«фильтр».
Можно добиться сколь угодно точного вычисления коэффициентов видимости при стремлении к нулю диаметра элементов разбиения ячейки-цели при численном интегрировании. Использовались итерации с пересчетом коэффициентов переноса для тех источников, для которых не были выполнены контрольные соотношения (93) с достаточной точностью. Увеличивалось число счетных лучей источник-цель, и вычисление плохо посчитанных коэффициентов проводилось вновь автоматически — до достижения необходимой точности.
Вычисленные значения коэффициентов нормировались. Использование нормированных значений коэффициентов обеспечивало явное балансное описание переноса с заданной точностью.
Решение задачи В. Г. Рогачева
Перейдем к изложению результатов счета задачи В. Г. Рогачева [18], имеющей аналитическое решение как в 2-мерном варианте, так и в 3-мерных вариантах.
Внутрь сферической полости через часть Sx граничной поверхности стационарно вводится излучение, оно полностью отражается (по Ламберту) от части ^ поверхности и свободно покидает полость через часть Sg поверхности.
Заданы диффузионные граничные условия, [7, стр. 124],
■Г = 0,25асТ4 Т 0,5g, (94)
а — постоянная Стефана. На каждой части поверхности сферы диффузионным условиям придается соответствующий смысл. Основой для решения задачи служит постоянство величины плотности одностороннего потока J на всей поверхности сферы. Решая систему трех уравнений, находим три значения температуры. При Sx = S3 и qx = Зас/8 после выхода на стационарное решение получаем Тх = (3/2)1/2, Т2 = (3/2)1/4, Т3 = (3/4)1/4
Для переноса внутри полой сферы S аналитически определяются угловые коэффициенты: р^ = Sj/S. В этом случае приближенный счет можно вести без дополнительного разбиения счетных областей вдоль образующей, вычисляя частичную сумму с шагом тг/2 по углу вращения. При счете задачи [18], например с 7 ячейками равной площади: по одной ячейке на Sx и 83, 5 ячеек на 82 — контрольные соотношения для коэффициентов ^ргз = 1 выполнялись с 6 выводящимися значащими цифрами.
Расчет переноса, выполненный с двойной точностью, дал стационарное решение так же с 6 точными значащими цифрами для ячеек на Si и S2, с 5 точными цифрами — на S3, шестая цифра отличалась от точного значения на единицу.
Для сферы, аппроксимированной поверхностями вращения 100-звенной ломаной линии со звеньями равной длины, вычисление коэффициентов доводилось до выполнения контрольных соотношений (93) с 6 значащими цифрами — итерациями с более мелким разбиением вдоль образующей и по углу вращения и с пересчетом коэффициентов для фиксированных сеточных источников. Точность результата счета переноса была та же.
Так же завершился расчет переноса внутри сферы и в 3-мерных вариантах постановки задачи. Задачу можно усложнять, произвольно разбрасывая по поверхности сферы области задания входного потока и области задания нулевого обратного одностороннего потока. Решение не меняется по величине, лишь бы для суммарных площадей этих областей выполнялось условие Si = S3.
Напомним, что с использованием интегрального уравнения в новой форме счет велся по явной, абсолютно балансной схеме — за счет нормирования коэффициентов видимости.
Задача о камере-обскуре
Рассмотрим решение задачи о камере-обскуре на основе изложенных алгоритмов. Пусть для простоты есть некоторое панно — квадрат плоскости, покрытый изображением. Изображение, как при печатании, — множество точек с разной степенью черноты. На некотором расстоянии есть две параллельных плоскости: плоскость преграды с отверстием и плоскость экрана. Пусть для определенности расстояние от плоскости с изображением до плоскости преграды равно расстоянию от преграды до плоскости экрана. Отверстие определяется как «фильтр». Угловые коэффициенты для каждой точки изображения определяются с помощью телесного угла, под которым отверстие видно из центра точки изображения. Перенос доли энергии излучения для точки изображения на экран порождает появление точки перевернутого изображения на экране.
Глава 7
СЧЕТ ЗАДАЧИ О ВЫХОДЕ ИМПУЛЬСА
ИЗЛУЧЕНИЯ ИЗ ЦЕНТРА ШАРА С УЧЕТОМИЗОТРОПНОГО РАССЕЯНИЯ
Шар равномерно заполнен веществом, изотропно рассеивающим излучение. Радиус шара равен двум эффективным пробегам рассеяния, ltot = к = Г Д = 2. Внешняя граница шара свободно пропускает излучение, приходящее на нее изнутри. Поток энергии извне отсутствует.
В начальный момент времени, t = О, на первом шаге по времени, пт = 1, в центральной внутренней ячейке шара, Ут=ь изотропно излучается импульс, величина излученной энергии принимается за единицу, S^=1 = 1. Других источников излучения нет.
После некоторого периода установления поток Fp^t) энергии, покидающей шар через внешнюю поверхность, должен описываться экспоненциальной функцией времени: Fr(P) = ехр{—At}. При радиусе шара в два пробега А = 0,4925, [8, 10].
Шаг по времени характеризуется длиной пути света ст, с = 1.
Счет переноса лучистой энергии велся с шагом ст = 0,01, ст = 0,005. Результат счета мало зависел от смены шага по времени и практически не зависел от перехода к описанию матриц коэффициентов переноса (не базовых коэффициентов) с одинарной точностью.
Важно выполнять соотношения:
ст kXl hy ltot'
При кж = b,hx = 0,4 рассматриваемая сеть S\^ имеет матричную структуру размерности 5x5 — для полусферы; кд = 10, kv = 25, ки = кд + kv + 1 = 36, ко = 2 х ку + кд = 60.
При кх = 10, hx = 0,2 сеть S1/4 имеет матричную структуру размерности 10 х 10 — для полусферы; кд = 20, ку = 100, ки = = кд + ку + 1 = 121, ко = 2 х ку + кд = 220, см. рис. 4.
Расчеты выполнялись с двойной точностью на PC Pentium (667 MHz).
Сеть счетных направлений переноса связывает друг с другом все ячейки. Каждая граничная ячейка-источник связана
Рис. 7.1. Матричная сетка для сферы, полученная проективным преобразованием квадратной сетки (kg = 20, кс = 100)
с каждым из элементов деления по углу граничных ячеек при описании переноса проникающей части энергии (по 3 коэффициента на направление) и с элементами внутренних ячеек — при переносе долей выведенной энергии (по 5 коэффициентов на направление).
Группа из 3 коэффициентов такова: (1 — сх^к)е^к, сх^ке^к, nijk. Группа из 5 коэффициентов: (1 - aimkyimb aimkeimk, nimkl (1 ^imk^Uimk, ^гтк^гтк-
Связи внутренних ячеек аналогичны.
Массивы коэффициентов при числе к^ разбиений геометрии на элементы по углу вращения насчитывают к^ х 16250 чисел — при счете на сети 5x5.
Массивы коэффициентов при числе к^ разбиений геометрии на элементы по углу вращения насчитывают к^ х 134 400 чисел — при счете на сети 10 х 10.
При численном интегрировании точность счета каждого коэффициента и время его счета определяется числом кх х к2 разбиений поверхности граничной ячейки-цели: на кх частей вдоль образующей и на к^ частей по углу вращения.
В расчете каждая граничная ячейка-цель разбивалась вдоль образующей на 5 частей и по углу вращения — на 100 частей.
Из центра источника в центры каждой из частей разбиения проводился луч, всего 500 лучей для счета одного сеточного взаимодействия источник-цель. Для каждого луча вычислялся путь и ослабление с весом этого луча для всех пересекаемых ячеек.
Для части коэффициентов точность счета проверяется на основе выполнения соотношений типа (3). Контрольные соотношения выполнялись с 4-6 точными значащими цифрами.
При счете серий задач [11-12] было установлено, что при однородно заданном значении пробега ltot = 100 на сети ячеек расчетной области с диаметром, равным 1, значения коэффициентов ослабления с 4 значащими цифрами совпадают с экспонентой ослабления, вычисленной с помощью среднего расстояния, Су « exp^—Pij/ltot)-
Все коэффициенты были пересчитаны при ls = 100 на сети ячеек задачи. Для всех матриц сравниваемые величины совпадали почти всегда с 4 значащими цифрами. Иногда последняя из 4 цифр отличалась на 1. Матрицы, вычисленные точнее, давали при сравнении совпадение 5 значащих цифр.
Для данной задачи отсутствуют граничные источники, в численном моделировании переноса участвуют лишь коэффициенты Ртр етр Ртр Ртгы етп' етты Ртр Ртр
Счет этих коэффициентов для сети 5x5 при кх = 5, к2 = 100 длился 25 секунд.
Счет коэффициентов для сети 10 х 10 при кх = 5, к2 = 100 длился 6 минут.
Счет 2000 шагов по времени моделирования переноса для сети 5x5 длился 23 секунды.
Счет 4000 шагов по времени для сети Юх 10 длился 3 минуты.
Искомый экспоненциальный показатель А имел разные значения при разном числе к^ разбиений ячеек по углу вращения, влиявшем на точность учета эффекта запаздывания, см. табл. 1, 2.
Таблица 7.1
Результаты счёта задачи о шаре-рассеивателе на сетке 5x5
к^
1
2
2
4
6
8
т
0,01
0,01
0,005
0,01
0,01
0,01
А
0,4560
0,4928
0,4927
0,4924
0,4895
0,4878
Таблица 7.2
Результаты счёта задачи о шаре-рассеивателе на сетке 10 х 10
2
2
4
6
10
т
0,01
0,005
0,01
0,01
0,01
А
0,49426
0,49426
0,49654
0,49488
0,49253
При счете на сети 8 х 8, ^ = 8, fei = 5, ^2 = Ю0, т = 0,01 был получен результат А = 0,49251. Время вычисления коэффициентов 1 минута, время выполнения 2000 шагов по времени — 2 минуты.
Последний результат (10 х 10, к^ = 10, кх = 5 , к2 = 100, т = 0,01, А = 0,49253), был проверен при к^ = 1000. Взаимодействие каждой ячейки-источника, граничной или внутренней, с каждой граничной ячейкой-целью рассчитывалось при рассмотрении 5000 лучей, идущих из центра источника в центры элементов разбиения поверхности цели. Счет коэффициентов занял 1 час. Был получен результат А = 0,49252. Изменилась лишь 6 значащая цифра. Пять цифр совпадают с опубликованным точным значением [8], см. также [10].
Глава 8
РЕШЕНИЕ ТЕСТОВОЙ ЗАДАЧИ О ВЫХОДЕИМПУЛЬСА ИЗЛУЧЕНИЯ ИЗ ЦЕНТРАУСЕЧЕННОГО ЦИЛИНДРА
Использовалась возникшая возможность выполнения сравнительного счета в составе комплекса программ САТУРН [26]. Одна и та же двумерная нестационарная задача решалась независимо:
по новой программе совместно с газодинамикой,
по основной программе комплекса совместно с газодинамикой.
Основная программа комплекса использует аппроксимацию интегро-дифференциального уравнения. Имеется большой опыт эффективного счета разнообразных сложных задач. Тестовая задача, о которой пойдет речь, была решена по основной программе за несколько минут на PC Pentium (667 MHz) — совместно с газодинамикой. Для новой программы решение данной задачи было важным тестом.
В центре цилиндра, усеченного с двух сторон перпендикулярными плоскостями, однократно — на первом шаге по времени — центральной внутренней ячейкой высвечивается импульс энергии. Других источников излучения не было. Геометрию этой задачи мы условно называли «консервной банкой». Радиус цилиндра был равен половине образующей.
Задача обладает дополнительной симметрией относительно центральной плоскости, делящей расчетную область на две части, левую и правую. Расчет велся с расстановкой точек-источников в правой части, а точек-целей — в обеих частях. Счет велся на сети 10 х 10 для правой половины при к\ = 5, /с2 = 100. Было 200 внутренних ячеек-целей и 40 граничных, разбитых на 10 частей по углу вращения — для аккуратного описания запаздывания, кх х к^ счетных лучей — для вычисления каждого коэффициента источник-цель.
Учитывалось поглощение и переизлучение вещества — в предположении локального равновесия этих процессов. На границе задавалось условие: половина приходящего одностороннего
потока отражается обратно внутрь области расчета по закону Ламберта.
Параметры вещества были совершенно случайны. Задача не имеет прикладного значения. Важен был сравнительный счет по двум программам — совместно с газодинамикой.
По основной программе задача решалась с крупным шагом по времени, общим с шагом решения уравнений газодинамики
По новой программе было выполнено 20 000 шагов, столько же — по газодинамике.
За это время радиус центрального осевого сечения и расстояние вдоль оси увеличиваются вдвое. Форма граничной поверхности приближается к сфере. Место стыка цилиндра и дисков — донышек «консервной банки» — продолжало выделяться выступом в одну граничную ячейку с каждой стороны выступа. Выпуклость границы была нарушена. Использовался лучевой алгоритм определения видимости.
Сравнение результатов легче всего было проводить по форме границ. Сравнивались координаты узлов ломаной линии, третьи цифры результата отличались не более чем на 1-2 единицы.
Обнаружилась необходимость совершенствования программы для вычисления точек пересечения луча с ячейками внутренней сети. Если луч пересекает грань ячейки с образующим отрезком, имеющим номер по, в массивах Ml (по), М2(по) сообщаются номера разделенных гранью внутренних ячеек. Для граничных отрезков один из номеров примыкающих ячеек равен нулю. Если луч пересекал одну из двух ячеек, примыкающих к грани, ясно в какую ячейку перейдет луч после этого.
Однозначность перехода исчезает, если луч выходит пересекая ребро ячейки. К ребру примыкают четыре ячейки. Лишь одну из них пересекает луч при своем продолжении. Не было предусмотрено описание этой топологической роли узлов сети ячеек. Не указывались номера примыкающих ячеек, четырех ячеек — при матричной структуре сети ячеек. Поиск соответствующей информации происходил с появлением диагностики «не вычислена точка выхода из ячейки». Это значило, что по ошибке поиск пересечения луча происходит не с той ячейкой.
Перебор всех соседних с ребром ячеек приводил к успешному счету, но с большими затратами времени. Можно просто выбрасывать из рассмотрения счетный луч, пересекший ребро или грань — близко с ребром, так что возникала упомянутая диагностика. Этот прием много лет использовался в расчетах по методу Монте-Карло. Основой является малая вероятность
пересечения с ребром, по сравнению с вероятностью пересечения с гранью.
При выбрасывании «трудного» луча счет коэффициентов занимал 15 минут работы PC Pentium (667 MHz). При поиске правильного продолжения пути луча — 2 часа. Коэффициенты пересчитывались через 1000 шагов по времени. Счет 1000 шагов занимал 5 мин. работы PC. Чаще чем через 1000 шагов пересчитывать коэффициенты не требовалось в связи с малым изменением телесного угла для основного источника.
Заключение
Новое интегральное уравнение дает аналитическое описание переноса излучения через область для фиксированного точечного источника. На его основе возможно явное численное описание переноса с контролируемой заданной точностью.
В приближении серой материи успешно решены задачи, для которых имеются проверенные численные результаты.
В двумерном приближении используется широкий класс сеток.
Первая счетная схема на основе нового уравнения переноса получена в предположении квазистационарности границ области и коэффициентов взаимодействия излучения и вещества, т. е. в предположении неподвижности границ и стационарности коэффициентов взаимодействия за время преодоления светом диаметра области. Первая счетная схема наиболее эффективна для описания переноса в оптически тонких областях.
Очевидна возможность развития первой модели для описания переноса в оптически толстых областях. При исчерпании влияния фиксированного источника вдоль конкретного счетного луча необходимо обрывать счет для этого луча. Для повышения экономичности такой схемы счета необходимо научиться работать с неполными матрицами за счет формирования дополнительных справочных матриц адресов соседства для заданной оптической толщины области влияния каждого источника.
На первых тестовых расчетах достигалась точность результата в 6 и более цифр. Лет сорок назад, открывая первую отраслевую математическую конференцию трех организаций, ныне — Института прикладной математики им. М. В. Келдыша (Москва), ВНИИЭФ (Саров) и ВНИИП (Снежинск), — акад. Ю. Б. Харитон отмечал отсутствие для большинства математических расчетов сообщения о допусках точности полученных результатов в условиях проведенного численного эксперимента. Для разностных сеточных методик идеалом представляется двойной пересчет задачи с удвоением каждый раз числа ячеек сети. На практике это трудно выполнимо и, как правило, свидетельствует о порядке точности аппроксимации разностной методики, а не о точности счета по ней. Разве при удвоении числа точек величина
дисбаланса автоматически уменьшается? Расчет, как правило, проводится с использованием всех вычислительных средств — для достижения точности сохранения баланса энергии лишь в десятые доли процента. Вот в практике использования методов Монте-Карло тратятся дополнительные вычислительные усилия на статистическую оценку полученного результата.
Допуски точности измерения для экспериментальной техники определяются не в ходе уникального опыта.
По мнению автора, в новом интегральном подходе возможно, параллельно с расчетом, выполнять в условиях расчета решение тестовых задач с известным ответом. Примером такого расчета может явиться описание переноса излучения для вещества и границ, нагретых до постоянной температуры — расчет изотермии, рассматривавшейся Кирхгофом. Может быть, на основе такого расчета можно научиться делать автоматическую коррекцию средств, используемых в основном расчете, например, — коррекцию вычисленных коэффициентов распределения.
При проведении параллельных вычислений расчет каждого источника можно выполнять отдельно от других источников, имея в памяти соответствующую информацию о коэффициентах переноса во все цели и копию массивов ег> и ew. После этого нужно только объединить суммированием энергию, вынесенную на границу первой группой запаздывания для каждой копии массива ещ и аналогично — энергию, вынесенную во внутренние ячейки первой группой запаздывания для каждой копии массива ew.
Можно, наоборот, параллельно вести расчет коэффициентов переноса в сторону одной или нескольких целей. При этом в памяти процессора необходимы соответствующие коэффициенты, сведения о состоянии всех источников и соответствующие сокращенные варианты массивов ev и ew для описания запаздывания. На управляющей ЭВМ необходим расчет и рассылка сумм коэффициентов распределения — для параллельного нормирования коэффициентов, сборка и рассылка данных о состоянии источников — для начала параллельного счета нового шага по времени.
С неизбежными потерями точности по времени возможна разработка алгоритмов с использованием ослабленного предположения о квазистационарности коэффициентов взаимодействия излучения и вещества — с описанием продвижения излучения на несколько шагов по времени, а не на всю глубину расчета. Придется иметь дело с таким сложным объектом, как шаровой слой остановленного импульса фиксированного источника, его размещением на множестве сеточных ячеек.
Нет принципиальных трудностей для счета задач в спектральном многогрупповом приближении.
В заключение отметим возможность использования нестационарного расчета переноса излучения не только для звезд и атомного оружия — для диапазона температур источников излучения в несколько тысяч градусов. Приведем пример обыденного нестационарного процесса замыкания электропроводки.
Список литературы
Бугер П. Оптический трактат о градации света / Пер. с фр. — М., 1950. С. 7-9.
Peierls R. Ц Proc. Cambr. Phil. 1939. V. 35. Р. 610.
Ахиезер А., Померанчук И. Некоторые вопросы теории ядра. - М.-Л.: ОГИЗ, 1948.
Унзольд А. Физика звёздных атмосфер. — ИЛ, 1949.
Carlson B.G., Bell G.I. Solution of the Transport Equation by the SH Method / Proc. 2nd Intern. Conf. Peaceful Uses Atomic Energy. 1958. Geneva. P. 2386.
Владимиров B.C. 11 Tp. Матем. ин-та АН СССР. 1961. №61.
Зельдович Я.Б., Райзер Ю.П. Физика ударных волн и высокотемпературных гидродинамических процессов. — М.: Физматгиз, 1963.
Fleck LA. The Calculation of Nonlinear Radiation Transport by a Monte Carlo Method // Meth. Computational Phys. 1963. V. 1. P. 43-65.
Беллман P., Калаба P. Квазилинеаризация и нелинейные краевые задачи. — М.: Мир, 1968.
Рихтмайер Р.Д., Мортон К. Разностные методы решения краевых задач. — М.: Мир, 1972.
Зигель Р., Хауэлл Дж. Теплообмен излучением. — М.: Мир, 1975.
Мельниченко А.С., Огибин В.Н. Применение метода Монте- Карло к решению спектральных задач лучистого теплообмена Ц Журнал вычисл. мат. и мат. физ. 1977. Т. 17, №4. С. 1068-1074.
Бабаев Ю.Н., Баженов С.В., Дементьев Ю.А. Об особенностях решения одного интегрального уравнения переноса излучения // Вопросы атомной науки и техники. Сер. Методики и программы численного решения задач математической физики. 1978. Вып. 1(1). С. 7-9.
Трощиев В.Е., Юдинцев В.Ф. Итерационный метод постоянных поправок для решения спектральных задач переноса излучения // Вопросы атомной науки и техники. Сер. Мето-
дики и программы численного решения задач математической физики. 1978. Вып. 2. С. 13-16.
Дементьев Ю.А., Машинин Р.Ф., Миронова В.Ф., Ченцов Н. Н. Приближённое вычисление углового распределения для нестационарного теплового излучения с поверхности вращения // Вопросы атомной науки и техники. Сер. Методики и программы численного решения задач математической физики. 1983. Вып. 1(12). С. 32-34.
Андреев Е.С., Козманов М.Ю. Принцип максимума и минимума для решения задачи о распространении излучения в неподвижной среде // Вопросы атомной науки и техники. Сер. Методики и программы численного решения задач математической физики. 1983. Вып. 1(12). С. 35-36.
Баженов С.В., Дмитриев НА., Дибирова Л.М. и др. Методика численного решения двумерных задач нестационарного переноса излучения, обладающих симметрией вращения // Вопросы атомной науки и техники. Сер. Методики и программы численного решения задач математической физики. 1983. Вып. 3(14). С. 33-35.
Дементьев Ю.А., Машинин Р.Ф., Нагорный В. И. и др. Методика численного решения трёхмерных нестационарных задач энергопереноса излучением через область вакуума Н Вопросы атомной науки и техники. Сер. Методики и программы численного решения задач математической физики. 1984. Вып. 2. С. 26-28.
\Ч Гольдин В.Я., Юдинцев В.Ф. О расчёте квазирегулярного и регулярного режима для уравнения переноса // Вопросы атомной науки и техники. Сер. Методики и программы численного решения задач математической физики. 1985. Вып. 2. С. 43-49.
Андреев Е. С., Козманов М.Ю. Методика СВЕТ для решения задач теплового излучения // Вопросы атомной науки и техники. Сер. Методики и программы численного решения задач математической физики. 1985. Вып. 2. С. 50-57.
Garetts Е., Rudy Т.Е., Hickman R.B. A numerical solution to the steady-state equation of transfer in azimuthally symmetric geometries Ц JQSRT. 1985. V. 34, № 5. P. 417-433.
Гагарин А.П. Бугера-Ламберта-Бэра закон, Физическая энциклопедия. Т. 1. — М., 1988. С. 232.
Плетенева Н.П., Шагалиев Р.М. Аппроксимация двумерного уравнения переноса на четырёхугольных многоугольных пространственных сетках по разностной схеме с рас-
ширенным шаблоном // Вопросы атомной науки и техники. Сер.Математическое моделирование физических процессов. 1989. Вып.З. С. 34-41.
Drake D.J. A view-factor method for solving time-dependent radiation transport problems involving fixed surfaces with intervening, participating media // J. of Comp. Phys. 1990. V. 87. P. 73-90.
Лазарев Д.Н. Ламберта закон // Физическая энциклопедия. — М., 1990. Т. 2. С. 618-619. Анфимов Н.А. Лучистый теплообмен. Там же.
Беляков И.М., Бутнева О.В., Бухарова В.Я. и др. Организация комплекса САТУРН // Вопросы атомной науки и техники. Сер.Математическое моделирование физических процессов. 1992. Вып.З. С. 49-51.
Кукушкин А. Б., Коган В. И. Перенос излучения, Физическая энциклопедия. Т. 3. — М., 1992. С. 565, 566.
Babaev Yu.N., Bazhenov S. И, Basin A.A. et al. VNIIEV techniques and codes for solution of non-stationary 2D and 3D problems of radiation energy transport in domains of complex shape using viev factors. Conversion applications problems. VANT, Ser. Mathematical Modeling of Physical Processes, I 4, 1995. P. 38.
Basin A.A., Vasina E.G., Vatutin V. V. et al. Simulation of Radiant Energy Transport in ICE Targets in Diffusion-Vacuum Approximations. GSI-97-08, Report. High energy density in matter, produced by heavy ion beams. P. 68.
Zagrafov V. G. Study into Geometric Factors of Radiation Energy Transfer, ibid. P. 69.
Vatutin V. V., Zhitnik A.K., Morenko A.I. et at. Application of Monte-Carlo Method to Radiation Transport Calculation in Fusion Targets. JCI-report, 1997. High energy density in matter, produced by heavy ion beams.
Дойников A.C. Фотометрия, Физическая энциклопедия. T. 5. - М., 1998. С. 353.
Дементьев Ю.А., Карповцев Е.А., Морозова Е.В. и др. Метод коэффициентов ослабления // Вопросы атомной науки и техники. Сер. Математическое моделирование физических процессов. 2001. Вып. 2. С. 28-36.
Дементьев Ю.А., Карповцев Е.А., Морозова Е.В. и др. Дополнительное интегральное уравнение переноса излучения // Вопросы атомной науки и техники. Сер. Математическое моделирование физических процессов. 2001. Вып. 4. С. 35-41.
Дементьев Ю.А., Карповцев Е.А., Морозова Е.В. и др. Дополнительное интегральное уравнение переноса излучения // Труды РФЯЦ-ВНИИЭФ. 2003. Т. 4. С. 18-27.
Дементьев Ю.А., Кириллов А. И. Об одной аппроксимации лагранжевых поверхностей с помощью эйлеровой сетки в пространстве // Вопросы атомной науки и техники. Сер. Методики и программы численного решения задач математической физики. 1987. Вып. 1. С. 54-60.
Анисина И.М., Дементьев Ю.А., Жмайло В. А. и др. Распространение излучения в области осесимметричной формы от сферического объемного источника // 4-я Нижегородская сессия молодых ученых. Математика и математическое моделирование. Саров: СарФТИ, 1999.
Анисина И.М., Нарожная И. А. Распространение излучения в области осесимметричной формы от сферического объемного источника // Межд. конф, молодых ученых и студентов «Актуальные проблемы современной науки». Самара, 2000.
Научное издание
ДЕМЕНТЬЕВ Юрий Алексеевич
РАСПРЕДЕЛЕНИЕ ЛУЧИСТОЙ ЭНЕРГИИ ТОЧЕЧНОГОИСТОЧНИКА. НОВАЯ ФОРМА ИНТЕГРАЛЬНОГОУРАВНЕНИЯ ПЕРЕНОСА ИЗЛУЧЕНИЯ
Оригинал-макет: В.В. Худяков
Подписано в печать 10.06.05. Формат 60x90/16.
Бумага офсетная. Печать офсетная. Усл. печ. л. 8.
Уч.-изд. л. 9,6. Тираж 400 экз. Заказ №
Издательская фирма «Физико-математическая литература»МАИК «Наука/Интерпериодика»
117997, Москва, ул. Профсоюзная, 90
E-mail: fizmat@maik.ru, fmlsale@maik.ru;
http://www.fml.ru
Отпечатано с готовых диапозитивовв ФГУП «Производственно-издательский комбинат ВИНИТИ»140010, г. Люберцы, Московская обл., Октябрьский пр-т, 403
УДК 535.311+519.642.2
ББК 22.193; 22.343 f
ДЗО
#-
Издание осуществлено при поддержке Российского фонда фундаментальных исследований по проекту 04-01-14058д
Дементьев Ю. А. Распределение лучистой энергии точечного источника. Новая форма интегрального уравнения переноса излучения. - М.: ФИЗМАТЛИТ, 2005. - 128 с. - ISBN 5-9221-0620-1.
Вместе с известным уравнением Пайерлса новое уравнение замыкает интегральное описание переноса излучения через ограниченную область пространства. Новое интегральное уравнение описывает распределение энергии фиксированного точечного источника излучения по направлениям расходящихся лучей, распределение внутри объема вещества выведенной части энергии, распределение по граничной поверхности проникающей части энергии.
Численное моделирование выполняется в классе естественной гладкости решения по полностью балансной явной схеме. Счет коэффициентов распределения квазианалитичен. Счетная погрешность при вычислении коэффициентов контролируется локально и корректируется с помощью нормировки также локально. Обеспечены условия для глубокого распараллеливания вычислений.
Изложен тридцатилетний опыт решения большого числа нестационарных трехмерных задач внутри полостей сложной формы и первый опыт счета двумерных тестовых задач с учетом взаимодействия излучения и вещества.
Книга печатается в авторской редакции.
ISBN 5-9221-0620-1
© ФИЗМАТЛИТ, 2005
© Ю. А. Дементьев, 2005
ОГЛАВЛЕНИЕ
Предисловие автора 5
Введение 10
Глава 1. Интегро-дифференциальное уравнение переноса излучения 47
Уравнение в частных производных 47
Характеристическое уравнение 50
Алгебраическое уравнение 51
Глава 2. Интегральное уравнение Пайерлса 61
Глава 3. Новая форма интегрального уравнения переноса излучения 66
Первая форма интегрального описания распределения лучистой энергии точечного источника 66
Вторая форма интегрального описания распределения лучистой энергии точечного источника 71
Обсуждение результатов 75
Глава 4. Численная модель нестационарного переноса излучения в двумерном приближении для широкого класса сеток 81
Сеть шагов по времени. Сеть граничных и внутренних ячеек 81
Основные счетные величины описания вещества и излучения на шаге по времени 82
Распределение импульса излучения сеточного источника . . . 84
Замыкающие вычисления 86
Алгоритм учета запаздывания 87
Базовые коэффициенты и коэффициенты переноса 89
Счет переноса энергии на одном шаге по времени 91
Обсуждение результатов 91
Глава 5. Вычисление коэффициентов ослабления 93
Коэффициенты — двукратные определенные интегралы . . . . 93
Алгоритм пересечения лучом внутренних ячеек 99
Лучевой алгоритм видимости 101
Глава 6. Численное решение трехмерной задачи энергопереноса излучением внутри полости сложной формы с учетом запаздывания и затенения 103
Аппроксимация граничных поверхностей. Граничные условия 103
Аппроксимация интегрального уравнения Пайерлса 105
Аппроксимация интегрального уравнения в новой форме . . . 108
Особенности решения трехмерных задач НО
Решение задачи В. Г. Рогачева 112
Задача о камере-обскуре 113
Глава 7. Счет задачи о выходе импульса излучения из центра шара с учетом изотропногорассеяния 114
Глава 8. Решение тестовой задачи о выходе импульса излучения из центра усеченногоцилиндра 118
Заключение 121
Список литературы 124