Применение аппарата обобщенных функций для построения приближенных решений задачи переноса излучения Текст научной статьи по специальности «Математика»

УДК 519.63

Ю. И. Скалько, М. А. Мендель

Московский физико-технический институт (государственный университет)

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

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

Ключевые слова: приближенное решение задачи переноса излучения, обобщенные функции.

1. Введение

1.1. Постановка задачи

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

Рассматривается одночастотная краевая задача для стационарного уравнения переноса излучения. Задачу решаем в выпуклой области G £ R3, которая ограничена кусочно-гладкой границей dG из масса С1. Интенсивность равновесного излучения при заданной температуре среды 1Р (х). Полагаем, что на границе области задана плотность потока излучения Iqq (О, х), направленного внутрь области. Как частный случай, Iqq (О, х) = 0. О — единичный вектор, указывающий направление потока излучения.

Тогда для интенсивности потока излучения I (О, х) в направлении О в заданной точке пространства с координатами х можем сформулировать краевую задачу:

Где п — внешняя нормаль к границе области дС.

В работе [1] показано, что эта задача имеет единственное решение в классе непрерывных кусочно-дифференцируемых по пространственным переменным функций I (О, х).

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

(1.1)

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

Другая группа алгоритмов основана на симметризации задачи (1.1). А именно, от задачи (1.1) мы переходим к задаче для новой промежуточной переменной V(О, х) с симметричным оператором [1], [5]. Решение исходной задачи (1.1) получается из дифференциальной связи I (О, х) = V — - (ОУ) V. Выполняя дискретизацию полученной задачи с симметричным оператором для переменной V (О, х), получаем СЛАУ с симметричной матрицей. Однако такую симметризацию можно выполнить, только если 0 < ко ^ к (х). Кроме того, необходимость вычислять I (О, х) посредством дифференциальной связи снижает качество алгоритма.

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

2. Построение численного алгоритма

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

2.1. Предварительные сведения. Обобщенные функции

Пространство основных функций И = И (Ка) — все финитные, бесконечно дифференцируемые в Кп функции. Например, функция-шапочка:

Се выбирается из уеловия / ш£ (х) йх = 1.

Пространство обобщенных функций V = V (Ка) — линейные непрерывные функционалы на пространстве основных функций И (Яп)- Например, регулярные обобщенные функции, функционалы, определяемые локально интегри руемыми в Кп функция ми / (х) по формуле (¡, <р) = / / (х) (х) йх, <р (х) е Б (Кп).

2.2. Свертка обобщенных функций

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

Пусть / (х) и д (х) — локально интегрируемые функции в Кп, причем функция ^ (х) = / 19 (у) / (х — У)1 ^у также локально интегрируема в Кп. Сверткой д * / этих функций называется функция

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

Свертка обобщенной функции с функцией шапочкой ше (х) называется регуляризацией обобщенной функции / и /£ = / * ше = (/ (у) , ше (х — у)) е С^ (Кп)-

3. Численный алгоритм

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

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

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

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

Определим обозначение, которое будем использовать в дальнейшем. Пусть некоторая функция д (х) терпит разрыв вдоль некоторой поверхности L и п — единичная нормаль к этой поверхности. Будем обозначать через [д] — скачок функции д (х) при переходе через поверхность разрыва L вдоль норм али п, то ест ь [д] = д- —д+. д- — значение фун кции д (х) с той стороны поверхности L, из которой «выходит» единичная нормаль п, д+ — значение функции д (х) с той стороны поверхности L, в которую «входит» единичная нормаль п.

Лемма. Пусть I (О, х) — решение краевой задачи (1.1). Тогда справедливы следующие утверждения:

1) функция (OV) 1может терпеть разрывы, только там, где терпят разрывы функции kl или klp, и (OV) I + kl — [klp] = 0;

2) всюду, где I (О, х) - гладкая по пространственным переменным функция, она удовлетворяет уравнению (—OV + к) (OVI + kl) = (—OV + к) (klp);

3) на границе dG облает,и G выполняются соотношения — (OV) I+klac = klp, х e dG, (On) < 0 (OV) I + kl = klp, x e dG, (On) ^ 0.

Доказательство. Первое утверждение леммы следует из равенства (OV) I + kl — klp = 0, которое выполняется то всех точках гладкости I (O, х).

Второе утверждение леммы следует из того же равенства после применения к обеим его частям оператора (—OV + к).

Третье утверждение следует из граничных условий

I = Iqg, х e dG, (On) < 0, (OV) I + kl — klp = 0, ж e dG.

Краевой задаче (1.1) поставим в соответствие уравнение для непрерывных обобщенных функций I (O, х) e D' (Rn), справедливое для произвольной пробной функции из пространства основных функций (р e D (G):

(QVI + kl — klp) • (QVp + к<р) dV — J J (Qn) (1 + k)(I — I9G) pdS = 0. (1.2)

G dG, (Пп)<0

Теорема. Пусть I (O, x) — решение краевой задачи для, системы уравнений (1.1). Тогда, I (O, х), рассматриваемая как регулярная обобщенная функция I (O, х) e D' (Rn), удовлетворяет уравнению (3). В классе непрерывных на G обобщенных функций такое решение уравнения (3) единственно.

Другими словами, уравнение (3) имеет единственное решение в классе непрерывных на С обобщенных фикций, и это решение совпадает с I (П, х).

Доказательство. Существование. Покажем, что функция I («, х), являющаяся решением задачи (1.1), удовлетворяет в смысле обобщенных функций уравнению (1.2).

Пусть ПУ/ терпит разрыв на поверхности Ь. Воспользовавшись формулой Грина интегрирования по частям, для произвольной пробной функции из пространства основных функций р £ И (С) имеем:

\

с

(«VI + кI - к1р) ■ («Ур + кр) <Ш 11 («п)(1 + к) (/ - 1дрйБ

дС, (Пп)<0

р ■ (-«V + к) («VI + кI - к!р) дУ + ^

+

V

с

+Иоя

ь

+

«V/

+

к!

- [к1р]) рйБ+ «п («VI + к1дс - к1р) рйБ+

дС, (0п)<0

+ 11 «п («VI + к! - к!р) рйБ

\ ас, (Пп)>о Воспользовавшись леммой, получаем

с

(«VI + кI - к!р) ■ («Ур + к<р) дьУ - { { («п) (1 + к) - ^с) фйв = 0.

ЭС, (Пп)<0

Единственность. Пусть существует непрерывная, кусочно-дифференцируемая функция I* («, х) = I («, х), удовлетворяющая в смысле обобщенных функций уравнению (1.2). Тогда и («, х) = I* («, х) - I («, х) = 0 удовлетворяет однородному уравнению / / / («Уи + ки)-(«Ур + кр) (IV- / / («п) (1 + к) ирйБ = 0 для произвольной проб-

С до, (Пп)<о

ной функции р £ И (С).

Воспользовавшись формулой Грина интегрирования по частям, получаем

р ■ (-«V + к) («Уи + ки) (IV +

с

+

! J «п ([«Уи] + [ки]) рйБ+ ь

+ f f «п («Уи - и) рйБ+

дС, (Пп)<0

+ J J «п («Уи + ки) р(1Б \ до, (Пп)>0 )

Ввиду произвольности пробной функции р £ И (С), последнее уравнение означает, что функция V = «Уи+ки всюду непрерывна, а там, где гладкая, она удовлетворяет уравнению и граничному условию

Г -«Уи + ки = 0, \ V = 0, х £ дС, («п) ^ 0.

0

Следовательно, V = 0, ж £ С. Последнее в свою очередь означает, что непрерывная функция и всюду, где гладкая, удовлетворяет уравнению и граничному условию

( ОУи + ки = 0,

\ ОУи — и = 0, ж £ дС, (Оп) < 0.

Или, что то же самое, уравнению и граничному условию

( ОУи + ки = 0, \ и = 0, ж £ дС, (Оп) < 0.

Тогда и = 0, х £ С, что противоречит сделанному предположению: и (О, ж) = I* (О, ж) — I (О, ж) = 0. Этим противоречием доказывается единственность решения уравнения (1.2) в классе непрерывных на С обобщенных функций.

О

Уравнение (1.2) будем решать численно. Для этого выполним триангуляцию области С, тем самым разобьем область С на тетраэдры О^, которые перенумеруем индексом £ = 1:1: Т. В каждом тетраэдре вершины перенумеруем индексом I = 1:1:4. Обозначим через хь1 координаты вер шины I в тетраэдре Также индексом р = 1:1: Р перенумеруем все узлы сетки, полученной в результате триангуляции, и обозначим через хр координаты узла р этой сетки. Проделанное позволяет построить 3-индексный массив Ур1 следующим образом: Ур1 = 1, если хь1 = жр, и Ур1 = 0 - иначе.

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

(0, ж £ Нп (ж) = \ Щг ) = #.

Здесь мы использовали обозначение Ь = 1 : Т, = 1 : Т, I = 1:4, I' = 1:4; 5, 5¡'

1

0

Это построение можно выполнить следующим образом. В координатном пространстве £ построим тетраэдр с вершинами Р1 = (0, 0, 0), Р2 = (1, 0, 0), Р3 = (0, 1, 0) Р4 = (0, 0, 1).

Построим четыре линейных полинома Н1 (£) из условия Н1 = ■

Для каждого тетраэдра построим линейное преобразование £ =И*х, переводящее координатное пространство х в координатное пространство £ из условия =Оьхь1. Тогда полиномы Щ1 (х) = Н1 (^х) и будут искомыми базисными интерполяционными полиномами с узлами интерполяции х*1.

В дальнейшем всюду при записи выражений с индексными объектами будем придерживаться соглашения о суммировании. А именно, если в индексном выражении индекс встречается 2 раза, причем, как правило, один вверху, другой внизу, то подразумевается свертка-суммирование по всему диапазону значений индекса. Если же свертка не подразумевается, то это будет оговариваться особо или рядом с этими индексами будет ставиться о. Например, в выражении Н¿0 1 иЬо 1 по индексу I предполагается свертка, по индексу £ свертки нет.

х

функций, заданных на иС^ в виде I (х, О) = Щ1 (х) Ур11р (О).

Введя обозначение Нр (х) = Щ1 (х) Ур* 1, можем записать I (х, О) = Нр (х) 1Р (О).

Возьмем в качестве проверочных функций из пространства основных функций И (С) функции ц?' (х) = Н^1' (х) У^ у; где Н^1' (х) — регуляризация базисных функций И1'1' (х) = 5г'г51'1 Нц (х), т.е. свертка с функцией шапочкой; а Ур = рУр1 ■

Напомним, что бр 1, 5Р р - символы Кронекера.

Подставляем их в уравнение (1.2). Если е взять достаточно малым, то приводимые ниже выражения будут выполняться с любой, наперед заданной точностью:

V р

vt' v

Gt

(OV Ht I + kHt i) (OVH*'l' + kHf/ ^ dV--On I I (1 + к) Нь'1'HtidS

VptlIp —

— vp = vt'l'

dG, (Пп)<0

( J J klp (toVH* l' + kH1'1') dV

Gt

-On J J Iqg (1 + к) H* l'dS

\ dG, (Пп)<0 /

Таким образом, для определения Ip (O) мы получили СЛАУ Ар (O) Ip = fp'.

Матрица этой СЛАУ симметрична и положительно определенная. Решив эту СЛАУ, находим Ip (O) и строим приближенное решение задачи (1.1) I (O, x) = Нц (x) Ург 1 Ip (O).

4. Тестовый пример

В качестве тестовохх) примера для исследования предложенного а.;п'оритма рассмотрим следующую задачу: в области с коэффициентом рассеяния kout — 1 и интенсивностью равновесного излучения, равной 0, нжодится излучающий шар. Радиус шара R — 0.3, коэффициент рассеяния внутри шара kin — 10, интенсивность равновесного излучения внутри шара 1Р — 10. Необходимо найти плотность интенсивности излучения в любой точке, лежащей вне или внутри шара. Эта задача взята в качестве тестового примера потому, что для нее может быть выписано аналитическое решение.

Построим такое аналитическое решение.

O

^ + kl — klp, OS

(1.3)

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

Рис. 1

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

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

Плотность интенсивности излучения задается формулой

(1.4)

и = Ц1

где dO = dtf sin ada — элемент телесного угла.

Так как I (O) = I (а) те зависит от то U = d<р f^ I (a) sin ada = 1 f^ I (a) sin ada.

Введем переменную £ = cos а, тогда U = 1 /_! I (£)d£. Пусть точка находится вне шара (рис. 1).

Нас интересует интенсивность излучения в точке А в направлении вА, составляющим

угол а с

ol. о!

г,

од

= К, г> К

Из (1.3) следует, что интенсивность излучения в точке С в направлении В А равна

1с = 1Р (1 - ехр [~кгп\вд|)) .

Поскольку интенсивность равновесного излучения вне шара равна 0, то интенсивность излучения в точке А равна 1а = 1с ехр(-кои^ — (С\).

Длины отрезков —О и В (С равны соответств енно: \в(( \ = 2^/ В2 — т2(1 — £2), I-С\ = — VК2 — т2(1 — С2)) следовательно:

1а = 1Р (1 — ехр (—кгп2^К2 — г2(1 — £2))) ехр (—коЫ — л/В,2 — г2(1 — £2))) . (1.5)

Формула (1.5) верна только для лучей, проходящих через шар, и точек А, лежащих вне шара. Поскольку через шар проходят только лучи, для которых ^ 1 — (у)2 ^ £ ^ 1, то

U(г) = -

/ Ч

Vi_(W

d£, г > R.

Рассмотрим теперь случай, когда точка А лежит внутри шара (рис. 2).

В { —Т ~ ■—■—i а \ i 1 'I

1 * V 0 1 1 >

Рис. 2

Из (1.3) следует, что для направления луча О = т^щ интенсивность в точке А равна 1а = 1р — ехр кгП\А1^\^ ^, а для противоположного направления луча О =

i

интенсивность = 1Р — ехр (—кгП1 А&. Легко проверить, что длины отрезков АЙ и

—& равны |А£^| = + К2 — т2(1 — (2), || = —г( + у/ГВ2—~г2(Т—~^2). В данном случае все лучи, приходящие в точку А, приносят ненулевую интенсивность излучения, поэтому, проинтегрировав 1д по всем углам, получим, что

1 Г1

и(Г) = 2]0 1р (1 — ехр (—к^г0 К2-гЦ1-2)) <%+

+2 /11Р (1 — ехр (кгпг0 е-^*2-^1-^) %,

следовательно:

и (г) = !11Р (1 — ск (кыгО е-к^п2-г2(1-2) ^ г < К.

Выпишем общее аналитическое решение для задачи об излучающем шаре:

и(г) = <

' 2 1 1,( 1 — е-к™(^ г > ^ 11р (1 — ск (кцпг{) г < К.

(1.6)

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

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

Рис. 3

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

На рис. 3 изображены изолинии плотности интенсивности излучения на плоскости, проходящей через центр излучающей сферы. Расчет проведен на сетке из 76 691 узла и соответственно для 6, 50, 110 и 302 точек квадратурной формулы [4| интегрирования но единичной сфере.

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

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

Т а б л и ц а 1

N 6 14 26 50 110 302 974

е 8,50Е-04 4,57Е-04 3,00Е-04 2,40Е-04 2,38Е-04 2,37Е-04 2,37Е-04

Рис. 4

На рис. 4 дан график погрешности, построенный в логарифмическом масштабе. Первые 4 точки ложатся на линию тренда, наклон которой равен 0,6. При числе узлов квадратурной формулы, равном 110, достигается предел сходимости на сетке из 76 691 узла. Таким образом, можно сделать заключение, что метод сходится по угловым переменным с порядком 0,6.

На рис. 5 приведена разность точного и приближенного решений на прямой, проходящей через центр сферы, при различном количестве узлов квадратурной формулы: егги6,егги2б,егги50,егги302 соответственно для 6, 26, 50, 302 и 974 узлов. Здесь точка 0 на оси абсцисс соответствует центру сферы, 0,3 границе излучающего шара, а 1 границе области расчета. Видно, что внутри шара погрешность достигает минимального значения уже при 6 узлах квадратурной формулы. В области большого градиента решения на границе шара имеем наибольшую погрешность, которая в большей степени обуславливается

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

Рис. 5

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

Из рис. 6 видно, что наибольшее отличие от точного решения наблюдается в области большого градиента функции, область [0,25; 0,35]. В остальных областях погрешность быстро стремится к 0 с увеличением количества узлов квадратурной формулы. Рассмотрим теперь сеточную сходимость метода. В табл. 2 приведены значения среднеквадратичной погрешности метода по всем узлам сетки е при различном числе вершин тетраэдров (узлов сетки) N и 110 узлах квадратурной формулы.

Рис. 6

Т а б л и ц а 2

N 2319 10706 17394 76691

е 59,20Е-4 16Д0Е-4 8,42Е-04 2,37Е-04

Рис. 7

На рис. 7 построен график погрешности в логарифмическом масштабе. Наклон прямой соответствует порядку сходимости метода и равен 0,93.

6. Заключение

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

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

Литература

1. Владимиров B.C. Математические задачи односкоростной теории переноса частиц /7 Труды МИЛИ СССР. 1961. С. 3 158.

2. Владимиров B.C. Уравнения математической физики. М.: Наука, 1981.

3. Галаиии М.П., Лукин В.В., Чечеттт В.М. Методы решения уравнения переноса излучения для астрофизических моделей: препринт / ИПМ им. М.В. Келдыша. М., 2010.

№ 59. С. 30.

4. Лебедев В.И. Квадратурная формула для сферы 131-го алгебраического порядка точности // Докл. РАН. - 1999. - Т. 366, Вып. 6. - С. 741-745.

5. Марчук Г.И., Агошков В.И. Введение в проекционно-сеточные методы. — М.: Наука, — 1981.

Поступила в редакцию 17.04-2013