авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 |   ...   | 10 | 11 || 13 | 14 |   ...   | 15 |

«Посвящается пионерам освоения космоса Предисловие Фактически в настоящее время закладываются основы решения фундамен- тальных проблем, связанных с ...»

-- [ Страница 12 ] --

Rb r Rt, Совокупность всех точек A(r) сферической оболочки образует открытую область G с нижней Gb и верхней Gt границами – сферическими поверх – ностями, имеющими радиусы Rb и Rt соответственно. Векторное поле всех направлений световых лучей s(A) в каждой точке A(r) образует множество = + – единичную сферу, где + и – множества (полусферы) – – 8.1. Математическая постановка задачи. Метод решения направлений s с [0, 1] и [1, 0], отвечающих восходящему и нисходящему потокам излучения соответственно. Отметим, что в сферической модели система отсчета зенитного угла отличается от плоской модели, где + отвечает нисходящему, а – восходящему потокам. Фазовый объем – задачи суть tot [G Gt Gb ] = {(r, s) : r [G Gt Gb ], s }.

Опишем подробнее только что введенные множества:

G = {r : Rb r Rt, 0 } ;

Gb = {r = rb : r = Rb, 0 } ;

Gt = {r = rt : r = Rt, 0 } ;

= {s : 0, 0 2} ;

+ = {s = s+ : 0 /2, 0 2} ;

= {s = s : /2, 0 2}.

Для удобства записи граничных условий вводим множества – фазовые – области b Gb + = {(r, s) : r = rb Gb, s + }, t Gt = {(r, s) : r = rt Gt, s }, метки которых выбраны для наглядности совпадающими с первыми буквами слов bottom (дно) и top (верх), а также множество граничных точек gr = Gt Gb.

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

Полную интенсивность монохроматического (при фиксированном ) или квазимонохроматического (при фиксированных и ) стационарного из лучения (r, s), где индекс —длина волны (ниже везде опускается), в точке A(r) с радиус-вектором r в направлении s находим как решение общей краевой задачи (ОКЗ) теории переноса K = F in, = F t, = R + F b (8.1) t b в фазовой области с линейными операторами: оператор переноса D (s, grad) + tot (r) = (8.2) + tot (r) ;

s 536 Глава 8. Сферические задачи с осевой симметрией интеграл столкновений – функция источника – B(r, s) S = sc (r) (r, s, s ) (r, s ) d s, (8.3) d s = d d ;

интегродифференциальный оператор K D S;

оператор отражения q(rb, s, s ) (rb, s ) d s, s +. (8.4) [R](rb, s) = Суммарная аэрозольно-молекулярная индикатриса рассеяния нормирована по условию (r, s, s ) d s = 2 (r, cos ) d cos = 1;

(8.5) tot (r) и sc (r) – пространственные распределения полного сечения взаи – модействия (ослабления, экстинкции) нейтрального (незаряженного) опти ческого излучения с веществом, заполняющим область G, и суммарного аэрозольно-молекулярного коэффициента рассеяния;

– угол рассеяния – – – определяется из соотношения cos = cos(s · s ) = cos cos + sin sin cos( ), если s = (, ), s = (, ). Запись a · b означает скалярное произведение векторов a и b.

Выражение функции источника (8.3) с учетом азимутальной симметрии решения ОКЗ (8.1) (r,,, ) = (r,,, 2 ) можно представить как B(r,,, ) = (r,, cos + ) + (r,, cos ) (r,,, ) d d, = sc (r, ) d где введены обозначения (при этом [0, ]).

cos + = + sin sin cos( + ), cos = + sin sin cos( ).

Функция F in (r, s) – плотность источников излучения, расположенных – внутри области G;

F b (r, s) и F t (r, s) – источники излучения на границах – Gb, Gt сферической оболочки, определенные для лучей s, направленных внутрь области G.

Оператор R описывает закон отражения излучения от подстилающей поверхности, расположенной на уровне нижней границы Gb ;

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

8.1. Математическая постановка задачи. Метод решения При R 0 (или = 0) имеем дело с первой краевой задачей (ПКЗ) теории переноса K0 = F in, 0 = F t, 0 = F b (8.6) t b для сферического слоя с неотражающими граничными условиями: либо с абсолютно «черными» поглощающими границами, либо с «вакуумными»

прозрачными границами.

Скалярная функция с векторными аргументами (r, s) = (A, s(A)) = (r,,, ) определяется как решение ОКЗ (8.1) или ПКЗ (8.6) в фазовой области G + Gt + Gb +, tot = t b, и образует скалярное поле. Изменение функции (r, s) в окрестности некоторой точки A(r) характеризуется вектором grad. Дифференцирование скалярного поля функции (r, s) от векторных аргументов равносильно определению производной (r, s) по направлению s в точке A(r). Производная по любому направлению s равна проекции grad на это направление, следовательно, = (s, grad ) = (s, ).

s Краевая задача (8.1) рассматривается при естественных, вытекающих из физики исследуемого процесса, ограничениях на коэффициенты, источники и граничные условия:

а) tot (r) и sc (r) – ограниченные, кусочно-непрерывные и кусочно – дифференцируемые функции в области G;

б) (r, s, s ) – непрерывная функция угла рассеяния = arccos(s · s ), – имеющая кусочно-непрерывные производные по каждой переменной;

в) оператор отражения R – равномерно ограниченный, т. е.

– q(rb, s, s ) d s 0 1;

г) рассеивающая и поглощающая среда внутри области G, а также среды на нижней Gb и верхней Gt границах области G – немультиплицирующие – (без размножения);

д) источники F in (r, s), F t (rt, s ), F b (rb, s+ ) – ограниченные, кусочно – непрерывные или финитные функции (например, внешний мононаправленный параллельный поток, узкий пучок или точечный источник описываются -функцией).

Отметим, что все изложение ведется для области G – сферической – оболочки. Задача для полной сферы сводится к задаче со сферическим слоем путем постановки граничных условий с отражением на границе Gb, имеющей сколь угодно малый радиус Rb. Эти условия описывают 538 Глава 8. Сферические задачи с осевой симметрией прохождение излучения через внутреннюю сферу, ограниченную Gb. Если внутри полость, то на Gb ставится условие «прострела». Если внутренняя сфера является рассеивающей или поглощающей средой, то на Gb вводится «условие отражения», учитывающее ослабление излучения внутри малой сферы. Такие средства расширяют прикладные возможности рассматриваемой модели переноса радиации, в частности, они позволяют включить ряд задач астрофизики и физики планет.

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

Фиксируем некоторое направление s. Проведем через точку A(r(0)) в этом направлении прямую так, что уравнение прямой можно записать в виде r() = r(0) s,, (8.7) где D (r()) – текущая точка на прямой, A(r(0)) – фиксированная точка на – – прямой, от которой отсчитывается сдвиг = |AD| вдоль прямой. С помощью таких прямых можно взаимно однозначно преобразовать точки области G в точки (A, ). Это преобразование переводит функции, измеримые в G, в функции, измеримые на вдоль прямых с направлениями s. В силу (8.7) можно записать (r s, s) (r, s) =.

s = Прямые линии (8.7) – пути, по которым движутся фотоны, – являются – – характеристиками линейного дифференциального оператора / s. Уравнение переноса теперь можно переписать так:

+ tot (r s)(r s, s) = B(r s, s) + F (r s, s) (8.8).

= 8.1. Математическая постановка задачи. Метод решения Уравнение (8.8) при известной правой части можно разрешить явно путем интегрирования по характеристике:

(r, s) = (r s, s) exp tot (r s) d + + B(r s, s) exp tot (r s) d d + 0 + F (r s, s) exp tot (r s) d d. (8.9) 0 Если задача с ненулевыми граничными условиями, т. е. хотя бы одна из функций ft или fb :

= ft, = fb t b отлична от нуля, то ее можно свести к задаче с нулевыми граничными условиями преобразованием вида 0 (r, s) = f (r s) exp tot (r s) d (8.10), а решение представить в виде суммы (r, s) = 0 (r, s) + d (r, s).

Функция 0 (r, s) отвечает прямому излучению от источника и будет иметь те же свойства, что и f (r, s), но несколько сглаженные благодаря экспо ненциальному множителю в (8.10). Для функции d (r, s), соответствующей многократно рассеянной, диффузной компоненте, получается задача со сво бодным членом в уравнении F1 (r, s) = F (r, s) + S0 (r, s). (8.11) Выделим рассеяние первой кратности 1 = D1 F1, (8.12) т. е. представим суммарное поле в виде суперпозиции (8.13) d (r, s) = 1 + d1.

В уравнении для d свободный член будет иметь вид Fd1 (r, s) = S1 = S D1 F1 (8.14) и в результате интегрирования вдоль луча s (действие оператора D 1 ) и по всем направлениям единичной сферы (действие оператора S) оказывается 540 Глава 8. Сферические задачи с осевой симметрией сглаженным (по сравнению с F, fb, ft ) и по пространственным и по угловым переменным. Разрывы в функции F1 (r, s) по угловым переменным приводят к разрывам Fd1 (r, s) по r. Но с увеличением номера кратности рассеяния эти разрывы сглаживаются. Однако, разрывы по r в коэффициентах tot (r), sc (r), (r, s, s ) проявляются в (r, s) и B(r, s) на всех итерациях. При разработке численного алгоритма уделяется специальное внимание локальным свойствам решения, что позволяет повысить точность решения и описания поведения решения в окрестностях особых точек.

Из формулы (8.9) интегрирования по характеристике следует, что диффе ренциальные свойства (r, s) определяются соответствующими свойствами функций B(r, s), tot (r), F (r, s) в пределах G и гладкостью границ, которая характеризуется дифференциальными свойствами (r, s). Пространственные и угловые производные функции (r, s) существуют и ограничены. В окрест ности касательных направлений s (к линиям или поверхностям разрыва коэффициентов tot (r) и sc (r)) производные имеют особенности вида 1,, r s |r r | |s s | если радиус кривизны в точках r конечен. На поверхностях разрыва с нулевой кривизной угловые производные имеют особенности вида ln |s s |, где s – касательное направление, а сама функция (r, s) имеет разрыв вдоль – этих направлений.

Функция (r, s) является абсолютно непрерывной вдоль лучей s и на множестве точек r, s, отличных от r, s, удовлетворяет условиям непрерывности по Гёльдеру:

1/2 1/ s) (r, s) A (r + r, s + r s.

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

Пусть имеет место представление U 0 = 0, = [u1, u ], u u+1, U =, u= и вдоль характеристики с направлением s в точке A(r) значение (r, s) = f (r s, s) exp tot (r s) d + + E(r s, s) exp tot (r s) d d, (8.15) 0 8.1. Математическая постановка задачи. Метод решения где через E(r, s) обозначена правая часть (8.8), определяется через значение f (r s, s) в граничной точке D(r s), отстоящей от A(r) на расстоянии. Разбивая интегралы по два с учетом того, что отрезок |AD| = |AB| + |BD| a 1, или [0, ] = [0, a] + [a, ], преобразуем выражение (8.15):

a (r, s) = f (r s, s) exp tot d + tot d + a a + E exp tot d d + E exp tot d d = a 0 0 a f (r s, s) exp tot d exp tot d = + a a + E exp tot d d + E exp tot d d a a 0 a exp tot d f (r s, s) exp tot d = + a a + E( ) exp tot d d exp tot d + a a a + E( ) exp tot d d.

0 Выражение в фигурных скобках можно привести к виду a (r as, s) = f (r s, s) exp tot (r ( + a)s) d + +a a E(r ( + a)s, s) exp tot (r ( + a)s) d d + 0 542 Глава 8. Сферические задачи с осевой симметрией и в итоге получим явную рекуррентную связь между значениями функции (r, s) в текущей и предыдущей точках траектории луча s:

a (r, s) = (r as, s) exp tot (r s) d + a + E(r s, s) exp tot (r s) d d, (8.16) 0 которое совпадает с формулой интегрирования по характеристике (8.15), только вместо граничного значения f (r s, s) в (8.16) входит значение (r as, s) в точке B(r as), предшествующей расчетной точке A(r) и отстоящей от A(r) на расстоянии 0 a. При переходе от формулы (8.15) к формуле (8.16) использовано свойство аддитивности экспонент.

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

8.1.2. Интегрирование уравнения переноса по характеристике без ин терполяции. В этих условиях интегрирование проводится по всей длине траектории луча s от расчетной точки A(r) до точки входа луча s в область G через верхнюю Gt или нижнюю Gb границы по формуле (8.15), где f и E источники на границе и внутри слоя. Реализация расчета (r, s) по формуле (8.15) при известных функциях E(r, s) и f (r, s) осуществляется по следующему алгоритму (общая схема).

1. Определение границы (Gb или Gt ), которую пересекает луч s, входя в слой G, и расстояние от точки A(r) до точки Q(r s) пересечения траекторией луча s этой границы.

2. Нахождение координат точки Q с радиус-вектором r() = r s = (r, ) и углов (, ), описывающих направление луча s в местной системе координат точки Q;

вычисление f (r s, s) = f (r,,, ) из соответствующих граничных условий.

3. Определение координат точки D с радиус-вектором r( ) = r s = (r, ) на расстоянии от точки A и углов,, описывающих направление s в местной системе координат точки D.

4. Нахождение значений коэффициентов tot (r s) в точке D(r, ) и расчет интенсивности источников E(r s, s) = E(r,,, ) в точке D(r, ) в направлении s = (, ) с координатами (, ) направления s в местной системе координат с центром в точке D.

8.1. Математическая постановка задачи. Метод решения 5. Вычисление интегралов tot (r s) d ;

(8.17) () = I() = E(r s, s) exp [ ( )] d (8.18) по квадратурным формулам с адаптивным выбором шага интегрирования, учитывающим структуру коэффициентов tot (r), источника E(r, s), положение расчетной точки A(r), направление луча s, масштаб длины отрезка траектории луча s, поскольку при реальных размерах сферической оболочки Земли (радиус 6370 км) значения варьируют в большом диапазоне: при высоте атмосферы, равной 100 км, 0 2300 км.

6. Расчет (r, s) по формуле (r, s) = f (r s, s) exp[ ()].

Такой алгоритм (естественно, модульный) позволяет выбирать направления траекторий лучей s независимо друг от друга, а пространственную разностную сеть произвольно. При заданных функциях E(r, s) (в достаточно свободной форме: явно в виде формул, таблично и т. п.) можно рассчитывать значения (r, s), например, для одного или нескольких направлений s в одной или нескольких точках A(r).

8.1.3. Интегрирование уравнения переноса по характеристике с интер поляцией. Сеточно-характеристический подход используется для расчета полного набора значений сеточных функций (rm, sn ) при известных сеточных функциях E(rm, sn ), например, в приближении многократного рассеяния, когда E(rm, sn ) является функцией источника и содержит сеточные значения интеграла столкновений. В этом случае в областях G и по всем перемен ным вводится разностная сеть. Радиусы rk и полярные углы l образуют пространственную разностную сеть, учитывающую геометрию задачи, т. е.

каждой точке слоя с радиус-вектором rm ставится в соответствие пара (rk, l ).

Совокупность возможных расчетных направлений луча sn в каждой точке слоя определяется парой углов (i, j ).

Определение интенсивности излучения (rm, sn ) в точке A(rm ) в направ лении sn проводится интегрированием уравнения (8.19) + tot (rm )(rm, sn ) = E(rm, sn ) sn 544 Глава 8. Сферические задачи с осевой симметрией по характеристике – лучу sn :

– (rm, sn ) = (r, sn ) exp tot (rm sn ) d + + E(rm sn, sn ) exp tot (rm sn ) d d. (8.20) 0 Сдвиг = |rm r | вдоль луча sn берется в пределах границ расчетной пространственной ячейки, представляющей собой кольцо, ограниченное двумя коническими (y1 и y2 ) и двумя сферическими (с радиусами r1 и r2 ) поясами.

Условно координаты точки A, в которой вычисляется (rm, sn ), обозначим через (r1, y1 ). Луч может войти в ячейку либо через границы r2 и r1 (лучи s2 и s4 ), либо через границы y2 и y1 (лучи s1 и s3 ). Оказывается, вдоль луча s4 на участках BC и AB величина имеет разные знаки;

= 0 в точке B.

Аналогично, на участках AD и DE луча s3 разные знаки имеет величина t.

Для построения более удобного алгоритма будем проводить интегрирование по лучу только на участках монотонного изменения переменных r и y, что соответствует постоянству знаков и t на этих участках. В таком случае в качестве точки с радиус-вектором r s (a = в формуле (8.16)) может быть получена одна из точек B, D, M, N. Обозначим координаты этой точки через (r, ), а направление s в этой точке – через (, ). Для вычисления – интеграла на участке характеристики [0, ] для функции E(r, s) вводится интерполяция по между значениями в узлах E(r,,, ) и E(r,,, ).

Значение (r,,, ) и E(r,,, ) при аргументах r,,,, отличных от узлов разностной сети, вычисляются с помощью интерполяции между соседними узлами разностной сети. Четыре возможных положения точки (r, ) порождает четыре случая интерполяции.

Реализация алгоритма интегрирования по характеристикам требует деталь ной проработки уравнений характеристик с тем, чтобы осуществлялась такая последовательность перебора аргументов rk, l, i, j, которая обеспечивает неперемешивание итераций при получении значений (rk, l, i, j ). Кроме того, необходимо уметь определять значения всех четырех переменных r,,, в любой точке траектории луча s по заданным значениям r1, 1, 1, 1. С помощью геометрических построений это сделать невозможно.

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

Алгоритм расчета одного луча sn = (i, j ) в точке сводится к следующей последовательности.

8.1. Математическая постановка задачи. Метод решения 1. Определение точки входа луча sn в расчетную ячейку по заданным значениям rk, l, i, j и координат r,,,.

2. Вычисление (r,,, ) и E(r,,, ) по формулам интерполяции между узлами разностной сети.

3. Расчет интегралов (см. формулу (8.20)) методом квадратур с адаптивным шагом, учитывающим, в частности, структуры задания коэффициентов tot (r, ). Отметим, что табличное задание tot (r, ) может быть более подробным, нежели пространственная разностная сеть, выбор которой зависит от оптической плотности. Если коэффициенты tot (r) являются ступенчатой функцией по и границы разрывности коэффициентов по входят в разностную сеть l, то лучше перейти к интегрированию по радиусу. Если имеются зоны с постоянными коэффициентами tot (r), то в пределах этих зон некоторые интегралы при реализации формулы (8.20) можно вычислить явно.

8.1.4. Основные компоненты решения. Общая краевая задача (8.1) ли нейная (по источникам) и ее решение можно искать в виде суперпозиции (аргументы (r, s) опускаем) = u + a + R = 0 + R, 0 = u + a.

Фоновое излучение 0 определяется как решение ПКЗ (8.6) и может содержать до трех фоновых компонент:

0 = t + b + in, 0 0 каждую из которых можно рассчитывать раздельно как решение ПКЗ с источниками F t, F b, F in соответственно:

Kt = 0, t = F t, t = 0;

(8.21) 0 0 t b Kb = 0, b b = F b;

= 0, (8.22) 0 0 t b Kin = F in, in in = 0, = 0. (8.23) 0 0 t b В каждой из этих компонент можно выделить по две составляющие:

t = t + t, b = b + b, in = in + in.

u a u a u a 0 0 Прямое излучение от источника на границе Gt – решение задачи – Dt = 0, t = F t, t = 0;

(8.24) u u u t b очевидно, что s = s, s = s+ +.

t t = 0 для для u u Отметим, что при освещении САЗ внешним мононаправленным потоком из-за непрозрачности Земли имеет место «область тени» Gten (ОТЗ), в 546 Глава 8. Сферические задачи с осевой симметрией которую не проникает прямое излучение солнечного потока, т. е. t = 0, u если (r, s) ten. Область тени Gten ограничивается сбоку цилиндрической поверхностью с радиусом Rb и осью OZ, сверху – вогнутой полусферической – поверхностью с радиусом Rb, а снизу – частью выпуклой сферической поверх – ности с радиусом Rt, которая попадает внутрь указанного выше цилиндра.

Через ten обозначаем фазовую область, включающую ОТЗ Gten и множества направлений s ten (r), зависящие от положения точки с радиус-вектором r Gten. «Область тени» существенно осложняет решение задачи (8.1) и рассматривается подробнее отдельно.

Прямое излучение от источника на границе Gb – решение задачи – Db = 0, b b = Fb ;

= 0, (8.25) u u u t b очевидно, что s = s+ +, s = s.

b b = 0 для для u u Прямое излучение от внутрених источников – решение ПКЗ – Din = F in, in in = 0, = 0. (8.26) u u u t b Суммарное прямое излучение от источников u = t + b + in u u u есть решение ПКЗ Du = F in, = Ft, = Fb. (8.27) u u t b В многократно рассеянном суммарном фоне – решении ПКЗ с нулевыми – граничными условиями = 0, =0 (8.28) Ka = Bu, a a t b можно выделить компоненты:

a = t + b + in, a a a соответствующие компонентам функции источников t b in Bu = Bu + Bu + Bu, Bu Su, Bu St, Bu Sb, Bu Sin, t b in u u u которые можно искать раздельно как решения следующего набора ПКЗ:

Kt = Bu, t t t = 0, = 0;

(8.29) a a a t b Kb = Bu, b b b = 0, = 0;

(8.30) a a a t b Kin = Bu, in in in = 0, = 0. (8.31) a a a t b 8.1. Математическая постановка задачи. Метод решения Решение каждой из задач (8.29)–(8.31) для диффузных компонент, которые можно записать в обобщенной форме как ПКЗ = 0, = 0;

(8.32) Kd = Bd, d d t b d = t, t для задачи Bd = Bu (8.29);

a b, b для задачи d = Bd = Bu (8.30);

a in, in для задачи d = Bd = Bu (8.31), a методом последовательных приближений – методом простых итераций по – кратности рассеяния есть сумма ряда Неймана:

N (8.33) d = n, n= слагаемые которого определяются из систем рекуррентных ПКЗ = 0, =0 (8.34) Dn = Bn1, n n t b с функцией источников Bn1 Sn1, B0 Bd. (8.35) Задача для определения подсветки R, обусловленной влиянием отража ющей подстилающей поверхности, – это ОКЗ – KR = 0, = 0, (8.36) R R = RR + ER, t b где источником инсоляции является освещенность (яркость, облученность) подложки, создаваемая фоновым излучением, ER (rb, s+ ) R0. (8.37) В зависимости от методических и прикладных интересов решение исходной ОКЗ (8.1) можно искать в виде суперпозиции = t + b + in (8.38) решений ОКЗ отдельно с каждым из источников инсоляции:

Kt = 0, t = F t, t = Rt ;

(8.39) t b Kb = 0, b b = Rb + F b ;

= 0, (8.40) t b Kin = F in, in in = Rin.

= 0, (8.41) t b 548 Глава 8. Сферические задачи с осевой симметрией При этом, как и в случае ОКЗ (8.1), в решениях ОКЗ (8.39)–(8.41), если необходимо, выделяем по три основные составляющие:

t = t + t + t = t + t, t = t + t ;

u a u a 0 R R b = b + b + b = b + b, b = b + b ;

u a u a 0 R R in = in + in + in = in + in, in = in + in.

u a u a 0 R R Фоновые вклады определяются как решения ПКЗ (8.21) для t, (8.22) для b, (8.23) для in. Прямое излучение – решение задача Коши (8.24) для – 0 t, задача Коши (8.25) для b, ПКЗ (8.26) для in. Диффузные вклады, u u u обусловленные многократным рассеянием в среде, – решения ПКЗ (8.29) для – t, (8.30) для b, (8.31) для in.

a a a Вклады подсветки от подстилающей поверхности находятся как решения следующих ОКЗ, аналогичных ОКЗ (8.36):

Kt = 0, t t = Rt + E t ;

= 0, (8.42) R R R R t b Kb = 0, b b = Rb + E b ;

= 0, (8.43) R R R R t b Kin = 0, in in = Rin + E in ;

= 0, (8.44) R R R R t b источники определяются подобно (8.37):

E t Rt, E b Rb, E in Rin.

0 0 Так что ОКЗ (8.36), (8.39)–(8.41) можно записать в обобщенной форме Kq = 0, = 0, (8.45) q q = Rq + E ;

t b q = t, E = Et для задачи (8.42);

R b Eb для задачи q =, E= (8.43);

R in, E in для задачи q = E= (8.44);

R для задачи q = R, E = ER (8.45).

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

Это обстоятельство позволяет включать в задачи (8.24)–(8.26) различные источники, не нарушая последовательности приближений.

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

8.1. Математическая постановка задачи. Метод решения Один из простых приемов заключается в том, что вся фазовая область G разбивается на подобласти Gnk = Gn k, в каждой из которых выбирается одна характерная фиксированная точка rn с направлением sk. По значениям интенсивности излучения в этой точке d1 (rn, sk ) проверяется выполнение критерия сходимости итераций в соответствующей подобласти. Для каждой подобласти Gnk вводится признак nk, по которому осуществляется выбор:

«заморозить» в Gnk результат предыдущей итерации или сосчитать новое приближение. Значение nk определяется при переходе на следующую итерацию либо условием для критерия сходимости, либо принудительно. Это позволяет на определенных итерациях исключать из рассмотрения отдельные подобласти Gnk, вклад в многократное рассеяние которых оказывается малым, или проводить расчет для полного фазового объема G (в частности, это необходимо для коррекции решения при завершении итераций).

Такой прием используется в методе редких сеток. Кроме того, для ускорения сходимости итерационного процесса применяется метод верхней релаксации с параметром в пределах от 1,5 до 2. Одним из приемов ускорения сходимости является «сглаживание» индикатрисы рассеяния.

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

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

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

8.1.5. К расчету однократного рассеяния. Если функция E(r, s) в фор муле (8.15) описывает в точке A(r) рассеяние прямого ослабленного излучения, дошедшего от источника, т. е. имеет вид (8.46) E(r, s) = F sc (r)(r, s, s0 ) exp [0 (r, s0 )], где 0 (r, s0 ) – оптическая толщина слоя между точкой A(r) и источником вдоль – направления распространения излучения источника s0, то кратко изложенный 550 Глава 8. Сферические задачи с осевой симметрией выше алгоритм интегрирования по характеристикам без интерполяции можно обобщить на случай селективного поглощения неселективного излучения, описываемого посредством функции пропускания. По существу, в этом алгоритме происходит суммирование вкладов в точке A(r) в направлении s от всех фотонов с направлениями s0, пересекающих s и испытывающих один акт рассеяния из направления s0 в направление s (в точке пересечения).

Если E(r, s) имеет место представление (8.46), то в интеграле I() (8.18) обе экспоненты можно объединить в одну и, выделив ослабление за счет рассеяния, записать I() = F sc (r s)(r s, s, s0 ) exp tot (r s) d exp [0 (r s, s0 )] d = = F sc (r s)(r s, s, s0 ) exp [sc ( )] P [w( )] d. (8.47) Процесс, описываемый интегралом I() (8.47), схематически можно изобразить так. Рассмотрим, например, случай внешнего параллельного потока в направлении s0. Фиксируем на направлении s точку C(r s), отстоящую на расстоянии от расчетной точки A(r). В точке C испытывает рассеяние фотон, пересекший границу Gt в точке B в направлении s0. В результате прохождения отрезков BC и CA происходит ослабление излучения за счет рассеяния, которое описывается функцией exp [sc( )], и за счет поглощения, которому соответствует функция пропускания P [w( )];

sc ( ) – суммарная – оптическая толщина отрезков BC и CA, определенная только по рассеянию;

P [w( )] – эмпирическая функция пропускания, зависящая от эффективного – содержания поглощаемой массы w( ) на пути фотона BC + CA. Величина sc (r s)(r s, s, s0 ) описывает увеличение интенсивности в направлении s в результате рассеяния луча s0 в точке C под углом = arccos(s · s0 ). Суммирование в интеграле I() осуществляется по всем точкам отрезка AQ таких, что 0.

Аналогично можно записать (8.48) exp [ ()] = exp [sc ()]P [w()], где оптическая толщина (по рассеянию) отрезка AQ sc () = sc (r s) d ;

(8.49) 8.1. Математическая постановка задачи. Метод решения эффективная масса поглощающего вещества на пути AQ w() = a (r s) d ;

(8.50) a (r) – эффективная концентрация поглощающей компоненты с учетом по – правок на зависимость от температуры и давления;

P [w()] – функция – пропускания, описывающая поглощение на отрезке AQ.

Если одновременно присутствует несколько поглощающих компонент (в том числе и поглощающие аэрозоли), то представляя суммарную функцию пропускания Psum в виде произведения (достаточно признанное приближение) (8.51) Psum = Pk [wk ] exp [a ], k можно вести Psum в формуле (8.46) и (8.48) и обобщить алгоритм на этот случай. В формуле (8.51) под знаком произведения стоят функции, описывающие поглощение через функции пропускания отдельных компонент Pk [wk ], а во втором множителе ослабление описывается коэффициентами поглощения, a – соответствующая оптическая толщина.

– Поскольку для () и I() (в общем случае) не удается получить результата в аналитическом виде, интегралы () и I() вычисляются приближенно методом квадратур с адаптивными сетками, учитывающими пространственную структуру оптических и метеорологических параметров, а также особенности функции E(r, s). Так границы разрывов параметров всегда включаются в узловые точки квадратурных формул. Особое внимание уделяется разрывам в источниках E(r, s), обусловленных границами тени (в области тени E(r, s) = 0, вне тени E(r, s) = 0), и -функциям, содержащимся в граничных функциях f (r, s).

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

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

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

8.1.6. К учету аэрозольного и молекулярного рассеяния. При модели ровании переноса излучения в земной атмосфере используются оптические параметры среды разной степени детализации. Кроме суммарных (обычно получаемых из натурных экспериментов) индикатрис рассеяния, особенно в ультрафиолетовой и миллиметровой области спектра длин волн, необходимо учитывать аэрозольные и молекулярные (рэлеевские) индикатрисы рассеяния с разной изменчивостью по высоте. Известно, что в земной атмосфере существует несколько аэрозольных слоев, которые различаются источниками их формирования, составом и оптическими свойствами. Аэрозольные индика трисы рассеяния либо измеряются в лабораторных условиях, либо рассчиты ваются по теории Ми, либо аппроксимируются функциями Хеньи-Гринстейна или другими функциями.

В качестве коэффициентов операторы содержат физические характери стики среды: пространственные распределения суммарного коэффициента ослабления (экстинкции) tot (r) = sc (r) + a (r), суммарных коэффициентов рассеяния sc (r) = a,s (r) + R (r) и поглощения a (r) = a,a (r) + R,a (r) с аэрозольными (a,s (r) и a,a (r)) и молекулярными (R (r) и R,a (r)) компонентами, суммарной индикатрисы рассеяния a,s (r) (r) a (r, ) + R (8.52) (r, ) = R (), sc (r) sc (r) которая определяется через индикатрисы аэрозольного a (r, ) и рэлеевского рассеяния (1 + cos2 );

ds = 2 (r, cos ) d cos = 1;

(8.53) R () = угол рассеяния из направления s = (, ) в направление s = (, ) находится из соотношения cos = + sin sin cos( ).

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

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

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

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

В аэрозольных индикатрисах рассеяния выделяем -анизотропию в виде:

a (r, ) = c(r)(1 cos ) + (r)a (r, ). (8.54) Полная исходная аэрозольная индикатриса рассеяния a, bx (r, ) перенорми руется по формуле a, bx (r, ) (8.55) a (r, ) = Aa, bx (r) = a, bx (r, ) d cos.

;

2Aa, bx (r) Как-то сглаженная (графически или путем аппроксимации функциями) исходная аэрозольная индикатриса рассеяния a, bx (r, ) также перенормиру ется:

a, bx (r, ) (8.56) a (r, ) = Aa, bx (r) = a, bx (r, ) d cos.

;

2 Aa, bx (r) Чтобы определить весовые коэффициенты в формуле (8.54), полную исходную аэрозольную индикатрису рассеяния записываем в виде a, bx (r, ) = f (r)(1 cos ) + a, bx (r, ) и ищем значение коэффициента f (r), исходя из условий нормировки (8.53), (8.55) и (8.56):

f (r) A (r) (1 cos ) + a, bx a (r, ). (8.57) a (r, ) = 2Aa, bx (r) Aa, bx (r) Откуда вытекает, что f (r) = Aa, bx (r) Aa, bx (r) 0, 554 Глава 8. Сферические задачи с осевой симметрией причем f (r) = 0, если «сглаженная» аэрозольная индикатриса рассеяния совпадает полностью с исходной аэрозольной индикатрисой рассеяния. Со поставляя выражения (8.54) и (8.57), получаем 1 (r) Aa, bx (r) f (r) (r) = c(r) =, =.

2Aa, bx (r) Aa, bx (r) Суммарная аэрозольно-молекулярная индикатриса рассеяния (8.52) с выделением -анизотропии представляется в виде, аналогичном (8.54):

(r, ) = v(r)(1 cos ) + (r, ) (8.58) с коэффициентом a,s (r) v(r) = c(r) a,s (r) + R (r) и сглаженной составляющей a,s (r) R (r) (r, ) = (r)a (r, ) + R ().

a,s (r) + R (r) a,s (r) + R (r) Если в суммарной индикатрисе рассеяния и суммарном коэффициенте рассеяния не выделять рэлеевскую компоненту, т. е. принять R (r) = 0, sc (r) = a,s (r), (r, ) = a (r, ), но оставить выделение -анизотропии, то в представлении (8.58) v(r) = c(r), (r, ) = (r)a (r, ).

Если анизотропию не выделять, то (r) = 1, c(r) = 0, (r, ) = a (r, ) = a (r, ).

Если выделить рэлеевское рассеяние, но -анизотропию не учитывать, то (r) = 1, c(r) = 0, v(r) = 0, (r, ) = (r, ).

Далее будем использовать наиболее общее представление суммарной аэрозольно-рэлеевской индикатрисы рассеяния a,s (r) (r) (r) c(r)(1 cos ) + a,s (r)a (r, ) + R (8.59) (r, ) = R (), sc (r) sc (r) sc (r) из которого можно получить частные случаи задания характеристик рассеяния среды.

С учетом представления суммарной индикатрисы рассеяния (8.58) и соотношения (1 cos ) = (s · s ) интеграл столкновений (8.3) можно условно записать в виде следующего выражения:

S S + Sc. (8.60) 8.1. Математическая постановка задачи. Метод решения Здесь введены обозначения для функции источника с сингулярной компо нентой индикатрисы рассеяния (8.58):

B S a,s (r)c(r) (r, s )(s · s ) d s = a,s (r)c(r)(r, s) (8.61) и со «сглаженной» индикатрисой Bc S c (8.62) (r, s ) (r, s, s ) d s.

Новую индикатрису рассеяния (r, s, s ) можно сделать достаточно гладкой и тем самым повысить точность квадратур при вычислении функции источника Bc (rm, sn ) (8.68). Аналогичный прием был использован Поттером в задачах астрофизики. Следует отметить, что на этом приеме основан метод подобия, когда от исходной задачи с оптической толщиной однородного слоя переходят к задаче с оптической толщиной, а весовой множитель c в (8.54) интерпретируется как эффективное поглощение, компенсирующее неполный учет индикатрисы рассеяния.

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

Рассмотрим задачу (8.21), когда источником излучения является внешний параллельный поток в направлении s0 = (0, 0 ) OZ, освещающий полусферу G0 Gt, r : r = Rt, G0 =, т. е. только часть G0 внешней поверхности Gt сферического слоя G. В этом случае F t = F (s0 ();

rb, s) = S (s · s0 ), где -функция определена по правилу 1, если s = s0, (s · s0 ) = 0, если s = s0 ;

функция Хевисайда («единичная ступенька») 1 для x 0, (x) = 0 для x 0.

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

0 =, 0 = cos.

s0 = (0 (), 0 ), 0 + =, 556 Глава 8. Сферические задачи с осевой симметрией Полагаем азимут 0 = 0 – азимут плоскости солнечного вертикала, про – ходящей через ось OZ и направление внешнего потока s0, т. е. азимуты направлений световых лучей s будем отсчитывать от плоскости солнечного вертикала.

Выделим прямое излучение u от источника, прошедшее через слой:

0 = u + a как решение задачи Коши (8.24). Многократно рассеянный фон a – решение – ПКЗ (8.29) с функцией источника Bu = Su.

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

Решение задачи Коши (8.24) интегрированием по характеристикам – лучам – s0 можно получить в явном виде:

(s · s0 )S exp [0 (r, s0 )], если r Gten, (8.63) u (s0 ();

r, s) = если r Gten.

0, Через Gten обозначена «область тени», внутри которой находятся точки r, координаты которых r, удовлетворяют неравенствам r cos (8.64) Rb, Rb r Rt,.

2 Граница BD описывается уравнением r cos (8.65) = Rb, из которого находится предельное значение зенитного угла lim, отвечающего точке D:

R R = b, lim = arccos b +.

cos lim (8.66) 2 Rt Rt Подобласть LDN M BDN M определяется независимыми для r и условиями Rb r Rt, lim.

В подобласти BDL BDN M при фиксированном значении = r радиусы точек K, Q, C изменяются в пределах Rb r R, R = Rb / sin r, а при фиксированном значении r зенитный угол для точек P, Q, R лежит в пределах R min = arccos b +.

min lim, r cos min = Rb, 2 r Оптическая толщина 0 (r, s0 ) вычисляется для «освещенных» точек A(r) Gten интегрированием вдоль траектории прямого луча s0 от расчетной точки 8.1. Математическая постановка задачи. Метод решения A(r) до предыдущей точки T (r s0 ) на луче s0, лежащей на пересечении лучом s0 верхней границы слоя Gt :

tot (r s0 ) d. (8.67) 0 (r, s0 ) = Для удобства изложения перепишем (8.63) в виде u (s0 ();

r, s) = (s · s0 )0, (s0 ;

r, s), (8.68) где если r Gten, S exp [0 (r, s0 )], 0, (s0 ;

r, s) = если r Gten.

0, В задаче (8.29) для многократно рассеянной компоненты a функция источника определяется через решение u (8.68) задачи (8.24) Bu = B0 Su = sc (r) (r, s, s )0, (s0 ;

r, s )(s · s0 ) d s = (8.69) = sc (r)(r, s, s0 )0, (s0 ;

r, s0 ).

Выделим вклад однократного рассеяния 1 :

N a = d = 1 + d1, где d1 = (8.70) n.

n= Этап 1. Полный вклад однократного рассеяния находится из ПКЗ = 0, =0 (8.71) D1 = B0, 1 t b с источником (8.69), содержащим полную аэрозольно-рэлеевскую индикатрису рассеяния. С учетом -анизотропии рассеяния (8.58) функция источника (8.69) преобразуется к виду:

B0 = B0, (s · s0 ) + B0, c, (8.72) где гладкая составляющая B0, c Sc u = sc (r) (r, s, s0 )0, (s0 ;

r, s0 ), (8.73) коэффициент в сингулярной компоненте B0, S u = sc (r)c(r)0, (s0 ;

r, s0 ). (8.74) Вследствие линейности краевой задачи (8.71) с учетом представления источника (8.72) можно разложить 1 на две компоненты:

1 = 1, (r, s0 )(s · s0 ) + 1, c. (8.75) 558 Глава 8. Сферические задачи с осевой симметрией Этап 2. Вклад однократного рассеяния со сглаженной индикатрисой рассеяния – гладкая составляющая фона – находится из задачи – – = 0, =0 (8.76) D1, c = B0, c, 1, c 1, c t b с гладкой функцией источника (8.73).

Этап 3. Сингулярная составляющая однократного рассеяния рассчиты вается с помощью решения задачи Коши = 0, (8.77) s = s0, D0 1, = B0,, 1, b где оператор D0 совпадает с оператором D при s = s0, т. е.

(8.78) + tot (r).

D0 = s В задаче для расчета d1 :

= 0, = 0, Dd1 = Sd1 + S1, d1 d t b проведем преобразование коэффициента экстинкции, переходя к представле нию (8.60) для функции источника:

= 0, = 0;

Dc d1 = Sc d1 + S1, d1 d t b оператор переноса Dc + [tot (r) a, s (r)c(r)]. (8.79) s Выделяем вклад второго рассеяния:

N d1 = 2 + d2, где d2 = n.

n= Этап 4. Вклад второго рассеяния находится как решение задачи = 0, =0 (8.80) Dc 2 = B1, 2 t b с функцией источника, которая определяется через полный вклад однократного рассеяния – решение задачи (8.71) по формуле (8.3):

– B1 S1 (8.81) или с учетом представлений (8.60) и (8.75) B1 = S 1 + Sc 1 = a, s c(r)1 + sc (r) (r, s, s0 )1, (r, s0 ) + Sc 1, c. (8.82) Учитывая -анизотропию в индикатрисе рассеяния (8.58) и однократном при ближении (8.75), в функции источника (8.81) можно выделить сингулярную компоненту:

B1 = B1, (s · s0 ) + B1, c. (8.83) 8.1. Математическая постановка задачи. Метод решения Гладкая составляющая функции источника B1, c a, s (r)c(r)1, c + sc (r) (r, s, s0 )1, (r, s0 ) + Sc 1, c ;

(8.84) коэффициент в сингулярной компоненте B1, a, s (r)c(r)1, (r, s0 ). (8.85) Исходя из представления источника (8.83), второе приближение можно записать в виде суперпозиции 2 = 2, (r, s0 )(s · s0 ) + 2, c. (8.86) Этап 5. Гладкая составляющая вклада второго рассеяния рассчитывается как решение задачи = 0, =0 (8.87) Dc 2, c = B1, c, 2, c 2, c t b с гладкой функцией источника (8.84).

Этап 6. Сингулярная компонента во втором рассеянии находится с помощью решения задачи Коши = 0, (8.88) s = s0 ;

D0, c 2, = B1,, 2, b оператор D0, c совпадает с оператором Dc (79) при s = s0 :

D0, c + [tot (r) a,s (r)c(r)]. (8.89) s Разложение (8.86) используется для расчета источника – сглаженного – интеграла столкновений B2, c Sc 2 = sc (r) (r, s, s0 )2, (r, s0 ) + Sc 2, c (8.90) в задаче, описывающей вклад диффузного излучения, начиная с третьей кратности рассеяния:

= 0, = 0. (8.91) Dc d2 = Sc d2 + Sc 2, d2 d t b Этап 7. Расчет вклада третьего рассеяния как решение задачи = 0, =0 (8.92) Dc 3 = B2,c, 3 t b с гладкой функцией источника (8.90).

Этап 8. Решение задачи (8.91) итерационно (N 4 – номер итерации) – по схеме:

Dc N = Sc N 1 + Sc 2, N N = 0, = 0. (8.93) d2 d2 d d2 t b В результате принятых мер по корректному учету сильной анизотропии в коэффициентах и источниках ПКЗ (8.21) для реализации этапов 1 необходима разработка четырех основных вычислительных модулей:

560 Глава 8. Сферические задачи с осевой симметрией I. интегрирование по характеристикам без интерполяции для расчета приближения однократного приближения;

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

III. интегрирование по характеристикам без интерполяции для расчета сингулярных компонент в направлении s0 ;

IV. расчет интегралов столкновений с гладкой индикатрисой рассеяния методом квадратур.

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

§ 8.2. Характеристики уравнения переноса Рассматривается перенос излучения в сферическом слое с осевой симметрией.

Пусть совокупность всех точек пространства A(r), r = (r, ), где r – – расстояние от центра сферы O, [0, ] – зенитный угол, отсчитываемый – от оси симметрии OZ, образует область G с границей gr, а совокупность возможных направлений распространения излучения s = (, ) – область – = [0, ] [0, 2] – единичную сферу = +, где полусфера – + = [0, /2] [0, 2] отвечает лучам s с [0, /2], или = cos [0, 1], полусфера = [/2, ][0, 2] – с [/2, ], [1, 0]. При азимутальной – симметрии () = (2 ) достаточно рассматривать интервал азимутов [0, ].

Значение интенсивности излучения находится как решение краевой задачи: (r, s) в фазовой области G удовлетворяют уравнению переноса (8.94) + tot (r)(r, s) = B(r, s) + F (r, s), s а на границе gr задаются краевые условия, зависящие от свойств границы и внешних источников;

B(r, s) – интеграл столкновений, F (r, s) – внутренние – – источники, tot (r) – коэффициент ослабления. В левой части (8.94) стоит – дифференциальный оператор переноса, содержащий производную по направ лению в точке A, которая для скалярной функции (r, s) = (A, s(A)) = (r,,, ) от векторных аргументов r, s совпадает с проекцией grad на это направление ((a, b) – скалярное произведение):

– = (s, grad ) = (s, = ).

s A s r 8.2. Характеристики уравнения переноса 8.2.1. Координатная запись дифференциальной части уравнения пере носа. Введем параметр – сдвиг вдоль луча – так, что – – (8.95) r() = r(0) + s, и положим по определению (r + s, s) =.

s r = Очевидно, что dr d d d = + + +.

=0 r d =0 d =0 d =0 d = Определив вид производных dr d d d,,,, d =0 d =0 d =0 d = мы тем самым найдем выражение для / s.

Из треугольника OAB, где OA = r(0), AB =, OB = r(), легко получить соотношения r 2 () = r 2 (0) + 2 2 r(0) cos( (0)), r() sin () = r(0) sin (0), c помощью которых тривиально находим dr d sin (0) = = cos (0) ;

.

d =0 d =0 r(0) Используя соотношения, выполняющиеся для траектории луча s:

cos = cos () cos () sin () sin () cos () = const, (8.96) r() cos () = r(0) cos (0) + cos, где – угол между лучом s и вертикалью z OZ, нетрудно найти – d sin (0) cos (0) =, d =0 r(0) d ctg (0) sin (0) sin (0) =.

d =0 r(0) В результате приходим к координатной записи sin ctg sin (8.97) = cos + cos.

s r r r Если вместо (8.96) взять выражение, которое получается при другой системе отсчета азимута:

(8.98) cos = cos cos + sin sin cos, 562 Глава 8. Сферические задачи с осевой симметрией то координатная запись производной примет следующий вид:

sin ctg sin (8.99) = cos cos +.

s r r r При нашей системе отсчета имеет место (8.96), поэтому мы будем использовать запись / s в виде (8.97). Нетрудно заметить, что (8.98) получается при сдвиге начала отсчета на.


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

dr rd rd rd = = (8.100) =.

cos sin cos sin ctg sin sin Из уравнения dr rd = cos sin получаем первый интеграл (8.101) r sin = c1 = const.

Два других члена rd rd = sin cos ctg sin sin дают второй интеграл (8.102) sin sin = c2 = const.

Из соотношения rd rd = sin cos sin с помощью (8.101) и (8.102) находим третий интеграл c+, если cos 0, c3 = c3, если cos 0, cos c+ = arcsin, (8.103) 3 c cos c = arcsin (8.104) +, 3 c а c4 введено для удобства:

1 c2 (всегда 0 1).

c4 = c При cos = 0 можно брать любое c3 (либо c+ либо c – однозначность 3– выбора вводится в алгоритме). Заметим, что c1 и c2 всегда неотрицательны, если,, принадлежат области, определенной выше;

очевидно c2 1.

Введем cos (8.105) cos x =, c 8.2. Характеристики уравнения переноса причем cos 0 (8.106) x = arccos.

c Тогда + arccos cos x = x (8.107) arcsin cos x =.

2 2 2 С помощью (8.107) легко определить область значений первых интегралов c+ и c при 0 :

3 c+ x x;

если 0, cos 2 c x x.

если 0, cos 2 Оказывается, что имеют место соотношения:

cos = c4 sin c+, если 0;

(8.108) cos если 0. (8.109) cos = c4 sin c3, cos Действительно, используя (8.103), получим cos = sin c+ cos + sin cos c+, (8.110) 3 c или (после элементарных выкладок) cos cos cos sin2 c+ 2 sin c+ + sin2 = 0.

3 c4 c Отсюда для sin c+ находим два выражения cos cos ± sin sin cos (sin c+ )1,2 =, 3 c одно из которых (со знаком «») дает (8.108), а другое отбрасываем, так как оно относится к области [, 2]. Аналогично получается и (8.109), если за исходное выражение взять (см. (8.104)) cos = sin c cos sin cos c. (8.111) 3 c 8.2.3. Геометрическая интерпретация первых интегралов оператора переноса. Выберем местную систему прямоугольных координат x, y, z с началом в точке (r, ) так, что ось z совпадает с z – линией, проходящей – через (r, ) и параллельной OZ, а оси x, y как-то ориентированы в плоскости Qhor, проходящей через (r, ) перпендикулярно к z. Введем единичные направляющие векторы n = {nx, ny, nz } для луча s;

n1 = {n1, n1, n1 } для xy z радиус-вектора r;

n2 = {n2, n2, n2 } для проекции r на плоскость Qhor (rhor );

xy z n3 = {n3, n3, n3 } для линии пересечения плоскости, соответствующей xy z = /2, с касательной плоскостью Qkac или что то же самое с плоскостью Qhor ;

n4 = {n4, n4, n4 } для проекции s на Qhor (shor );

n5 = {n5, n5, n5 } xy z xy z 564 Глава 8. Сферические задачи с осевой симметрией для перпендикуляра d к лучу s, проведенному из центра сферы т. O;

n6 = {n6, n6, n6 } для оси z.

xy z Легко видеть, что c1 = OC – кратчайшее расстояние s от центра. Если – ввести длину отрезка l = AC, то уравнение луча s можно записать так r = c1 n5 + l n, причем n n 5, так как s d. С другой стороны, как уравнение радиус-вектора, r = r n1. Следовательно, c1 n 5 + l n c1 n 5 + l n n1 = =.

r c2 + l В плоскости Qhor угол между линией и rhor равен /2, а потому n3 = n2, n3 = n2, n3 = n2 = 0.

x y y x z z Учитывая, что n n1 2 2 2 2 = y, x (n2 ) + (n2 ) = 1, (n1 ) + (n1 ) + (n1 ) = 1, x y x y z n2 n x y c1 n5 + lny c1 n5 + lnx y x n1 = n1 =,, x y r r 1 (n1 )2 = sin n1 = cos, z z (корень берется со знаком плюс, так как [0, ]), можно получить следующие выражения: 1 ny c1 ny + lny n2 = =, y r sin 1 (n1 ) z n1 c1 n5 + lnx x x n2 = =, x r sin 1 (n1 ) z c1 n5 + lny c1 n5 + lnx y x n3 = n3 =,.

x y r sin r sin Запишем двумя способами проекцию n на (n ):

1) сначала n проектируем на Qkac, а затем в касательной плоскости на – это дает – |n | = cos cos = sin sin ;

2 2) сначала n проектируем на Qhor, а затем в горизонтальной плоскости на – это дает – |n | = cos cos = sin cos, где – угол между shor и. Так как проекция n на при любом способе – проектирования должна быть одна и та же, то должно иметь место равенство sin cos (8.112) sin sin = sin cos, sin =.

sin 8.2. Характеристики уравнения переноса Далее определим через векторное произведение [a, b] [c1 n5 + l n, n] c sin = |[n1, n]| = = 1, (8.113) r r так как [0, ], cos = (n3 · n4 ) = n3 n4 + n3 n4. (8.114) xx yy Очевидно, что nx ny n4 = n4 = n4 = 0,,, x y z )2 )2 )2 ) (nx + (ny (nx + (ny (nx ) + (ny )2 = 1 (nz )2, 2 2 2 (nx ) + (ny ) + (nz ) = 1, 1 (nz )2 = sin, nz = cos, так как [0, ]. Подставляя значения проекций векторов n3 и n4 в (8.114), найдем c1 (nx ny ny nx ) 5 (8.115) cos =.

r sin sin Наконец, подставляя (8.113) и (8.115) в (8.112), получим sin sin = n5 ny n5 nx = [n5, n]z.

x y Поскольку n5 n, а n5, n – единичные векторы, то – [n5, n]z = cos, где – угол с осью z (или z) и линией T, являющейся векторным – произведением n5 и n и образующей с ними правую тройку векторов.

Принимая во внимание тот факт, что, [0, ], всегда имеем 0, sin sin = c а, следовательно, и cos 0. Тогда главное значение [0, /2].

Назовем плоскостью Q плоскость, проходящую через s, d и имеющую в качестве нормали линию T. Очевидно, что в этой плоскости лежит r и что она проходит через центр сферы O. В сечении сферы плоскостью Q получаем круг с радиусом сферы. Оказывается, что эта плоскость является координатной поверхностью, соответствующей азимуту луча s в точке A(r, ). Вращая s в точке A(r, ) по поверхности Q, получим набор всех [0, ] при фиксированных r,,. Введем = (/2 ) – угол наклона – оси z (равно как и z) к плоскости Q:

sin = sin = cos = sin sin.

Итак, интеграл (8.102) отражает очевидный теперь факт сохранения угла наклона оси z к плоскости Q, в которой лежит луч s.

566 Глава 8. Сферические задачи с осевой симметрией Помня, что так как [0, /2], то и [0, /2], можно определить 1 c2, cos = но это есть не что иное как c4, т. е. c4 = cos.

Нетрудно установить, что в касательной плоскости лежит линия T1, параллельная T и соответствующая азимуту + /2.

Действительно, плоскость Q проходит через r. T Q, существует T1 T и T1 r (T1 проводится через точку A(r, )). А такое T1 обязано лежать в Qkac, так как r Qkac и Qkac проходит через A(r, ). В силу того, что T s, так же T1 s. И если skac – проекция s на Qkac, то T1 skac, причем – skac ближе к линии = 0. Таким образом азимут T1 будет равен + /2.

Оказывается, что через соотношение (8.105) определяется косинус угла x между r и znp – проекцией оси z на плоскость Q.

– Действительно, спроектируем z на r двумя способами:

1) сначала на znp, затем на r – это дает n6 = c4 cos x;

– r 2) сразу на r – это дает n6 = cos.

– r Из условия равенства этих проекций получаем (8.105). Запишем область значений x :

0 x, если cos 0;

если cos = 0;

x=, x, если cos 0.

Как следствие соотношения (8.105), можно отметить, что а) при c4 = 1 (это возможно, если либо sin = 0, либо sin = 0, либо sin = 0 и sin = 0 одновременно) znp совпадает с z (x = );

б) при c4 = | cos | (когда cos = 0) znp или совпадает с r (x = 0), если cos 0, или направлена противоположно r (x = ), если cos 0;

в) при c4 = 0 (когда sin = 1, sin = 1, т. е. для = /2 и = /2 одновременно) znp обращается в точку, x – неопределенно;

этот случай – соответствует горизонтальной плоскости, проходящей через центр сферы O перпендикулярно оси z;

г) при c4 = | cos | (когда cos = 0) znp совпадает с линией азимута, если cos 0, или противоположно направлена, если cos 0, в Qkac, которая в этой ситуации параллельна оси z;

так что угол между znp и r при любом равен /2.

Для интеграла c3 = const можно дать две удобные геометрические интерпретации (обе в плоскости Q).

Первая интерпретация: c3 – угол между d и znp в Q-плоскости.

– Вторая интерпретация: c3 – угол, составляемый s c линией T, соответству – ющей линии пересечения плоскости Q с конусом = /2, вырождающимся в плоскость.

8.2. Характеристики уравнения переноса И та и другая интерпретация легко проверяются путем выражения угла 0 между d и znp через и x и сравнения c выражениями (8.103) и (8.104).

Отметим некоторые возможные характерные значения c3 и взаимное расположение s и znp, z при этих c3 :

cos c+ = c+ = 0, s znp, sz, sin c+ = 0 – возможно при x + = /2;

– cos c+ = 0 c+ = s znp 3,, sin c+ = 1 3 – возможно только при = x = 0;

– cos c+ = 0 c+ = s znp 3,, sin c+ = 1 3 – возможно при x + = ;

– cos c+ = c+ =, s znp, sz 3, sin c+ = 0 – возможно при x + = 3/2;

– cos c = c = 0, s znp 3, sin c = 0 – возможно при x = /2;

– cos c = 0 c = s znp 3,, sin c = 1 3 – возможно при = x;

– cos c = 0 /2, s znp – возможно при x = 0, =, – c =, sin c = 1 3/2, s znp – возможно при x =, = 0;

3 – cos c = c =, s znp, sin c = 0 – возможно при x = /2.

– Подчеркнем, что при sin c+ = 0 или sin c = 0 имеем s znp и s z, а 3 потому s Qhor.

c Что такое ?

cos c Легко установить, что это есть отрезок вдоль znp от центра сферы O до точки пересечения s и znp. Когда cos c3 = 0, s znp или s znp ;

568 Глава 8. Сферические задачи с осевой симметрией c при этом в и отражает тот факт, что s и znp обращение величины cos c не пересекаются.

Если cos c+ 0 или cos c 0, то точка пересечения s с znp лежит 3 на положительном направлении znp – этой точке на луче s соответствует – cos = c4.

Если cos c+ 0 или cos c 0, то точка пересечения s с znp лежит на 3 отрицательном направлении znp и этой точке соответствует cos = c4.

c дает отрезок вдоль линии T от центра сферы O до Величина sin c точки пересечения s с T. При sin c3 = 0 эта величина обращается в и, действительно, в этом случае нет пересечения s с T, так как s Qhor.

Если sin c+ 0 или sin c 0, то точка пересечения лежит на 3 положительном направлении линии T, а если sin c+ 0 или sin c 0, 3 то такая точка лежит на отрицательном направлении линии T.

c По существу, с помощью величины мы получаем точку на луче s, sin c соответствующую cos = 0.


В дальнейшем это подтвердится аналитическим исследованием уравнения для луча s.

8.2.4. Уравнения луча s в плоскости r, y (y = cos). Возьмем на тра ектории луча s две точки: фиксированную точку A[r(0)] с координатами r(0) = (r1, 1 ) и текущую точку D[r()] с координатами r() = (r, ), расстояние между которыми 0. Согласно (8.95) r cos = r1 cos 1 cos. (8.116) Пусть в точке A направление луча s определяется углами 1, 1, а в точке D –,.

– Рассмотрим случай cos 1 0 (c3 = c+ = const).

Умножим обе части (8.110) на r:

r cos = c4 cos c+ (r sin ) + c4 sin c+ (r cos ) 3 и примем во внимание, что r cos = ± r 2 c2.

r sin = c1, Тогда получим («+» при cos 0, « » при cos 0) r cos = c1 c4 cos c+ ± c4 sin c+ r 2 c2. (8.117) 3 Приравнивая правые части (8.116) и (8.117), после тривиальных выкладок получим одно из двух параметрических уравнений луча r 2 = 2 2b + c, (8.118) 8.2. Характеристики уравнения переноса где b = r1 cos 1, c = r1. При выводе использованы тождественные соотношения (8.108), (8.110).

Кривая r() является гиперболой. Действительно, её инварианты 0 b 0 1 = b2 c = c2 0, = b 0 c 1 = 1 0.

= 0 Равенство b2 c = c2 (8.119) устанавливается с помощью (8.110). Разрешим (8.118) относительно 1,2 = r1 cos 1 ± r 2 c2. (8.120) Какое значение r, cos и брать при фиксированных r1, 1 ?

Пусть r2 r1 r3. Тогда нетрудно установить, что если cos 1 0, c1 r2 (луч s1 ), то r2 c r = r2, r = r1 cos 1 r2 c2, 2 0;

cos = 1 r если cos 1 0, c1 = r2 (луч s2 ), то cos = 0;

r = r2, r = r1 cos 1, если cos 1 0, r2 c1 r1 (луч s3 ), то cos = cos 1 0;

r = r1, r = 2r1 cos 1, если cos 1 = 0 (луч s4 ), то r3 c r3 c2, cos = c1 = r1, r = r3, r = ;

1 r если cos 1 0 (луч s5 ), то r3 c r3 c2, cos = r = r3, r = r1 cos 1 + ;

1 r если cos 1 = 1, то r = r1 r2, c1 = 0, cos = 1;

r = r2, если cos 1 = 1, то r = r3 r1, cos = 1.

c1 = 0, r = r3, 570 Глава 8. Сферические задачи с осевой симметрией Второе параметрическое уравнение для луча r1 cos 1 cos (8.121) y() = 2 2b + c получается из (8.116) заменой r r() с помощью (8.118).

Уравнение для (y) можно привести к виду:

A 2 2B + C = 0 (8.122) с коэффициентами A = y 2 c2 sin2 c+, 4 B = r1 cos 1 y 2 r1 cos 1 c4 sin c+, C = r1 (y 2 cos2 1 ).

Детерминант этого уравнения det = B 2 AC = c2 y 2 (c2 y 2 ) 0, 1 в чем легко убедиться, если воспользоваться (8.111).

Разрешим (8.122) относительно :

если y 2 = cos2, то cos 1 y 2 cos 1 cos ± |y sin 1 | c2 y B ± det (8.123) y = = r1 ;

y 2 cos A если y 2 = cos2, то r1 (cos2 cos2 1 ) C y = = = 2 cos (r1 cos 1 cos r1 cos 1 ) 2B r1 (cos2 1 cos2 ) (8.124) =.

2c4 cos sin 1 cos c+ Вопрос о выборе знака в выражении для y, о выборе y и знака для c2 y sin2 c2 cos = ± =± (8.125) 1 y sin (для sin = 0, = 0 или, cos неопределен) гораздо сложнее, чем в случае r, и будет рассмотрен ниже в связи с классификацией различных возможных случаев прохождения луча s через область r, y.

Исключая из (8.118) и (8.121), придем к уравнению луча в общем виде:

r 2 y 2 2ry r 2 + = 0 (8.126) с коэффициентами c2 c = 1 24, 0, = cos2 cos r (cos 1 cos 1 cos ).

= cos 8.2. Характеристики уравнения переноса Это есть кривая четвертого порядка, у которой нас интересует лишь одна ветвь, соответствующая r 0.

Из (8.126) можно получить уравнение луча в виде:

если y 2 = cos2, то c1 c4 cos c+ y ± c1 cos c2 y r(y) =, y 2 cos причем если sin c+ 0, « + », если sin c+ 0, « », или по-другому y cos c+ | sin c+ | c2 y 3 (8.127) r(y) = c1 c4 ;

y 2 cos если y 2 = cos2, то c1 c r(y) = 2y cos c+ или в виде c1 c4 cos c+ ± cos r 2 c (8.128) y(r) =.

r Исследование уравнения луча s удобнее провести в форме (8.128).

Область значений r и y вдоль луча:

y2 c2.

0 c1 r +, При r = c1 cos c+ 0, 0, если + cos c+ = 0, y(c1 ) y0 = c4 cos c3 = 0, если (8.129) cos c+ 0.

0, если При r кривая y(r) имеет асимптоты c1 c4 cos c+ c ± cos 1 = ± cos.

lim y(r) = lim r r r r Значение y = 0 принимается при c r = r0 = (8.130).

| sin c+ | Действительно, должно быть c1 c4 cos c+ ± cos (r 0 )2 c2 = или c1 c4 cos c+ + c2 = (r 0 )2, cos 572 Глава 8. Сферические задачи с осевой симметрией отсюда c (r 0 )2 =.

sin2 c+ Но так как у нас r 0, то получаем как раз (8.130). В случае sin c+ = кривая y(r) нуль не пересекает.

Вычислим производную c1 c4 cos c+ r2 c2 ± c2 cos dy =.

dr r 2 c r2 dy обращается в – это точка Как видим, в точке r = c1 производная – dr dr поворота кривой y(r) или точка min r(y), поскольку здесь = 0.

dy Найдем точки rэ возможных экстремумов y(r) из условия dy = 0.

dr В случае знака «+» должно быть c2 c4 sin c+ c1 c4 cos c+ rэ c2 = 1 3 или c1 sin c+ = rэ c2 cos c+ 3 – это возможно, если sin c+ и cos c+ одинаковых знаков.

– 3 В случае знака «» должно быть c2 c4 sin c+ + c1 c4 cos c+ rэ c2 = 1 3 или c1 sin c+ = cos c+ rэ c 3 – это возможно, если sin c+ и cos c+ разных знаков. В обоих случаях – 3 c2 c 2 или rэ = rэ =.

| cos c+ | cos2 c+ При cos c+ = 0 y(r) экстремальные значения принимает при r.

Если y(r) берется со знаком «+», то при cos c+ 0:

если sin c+ 0, то 0 при r rэ, y(r) возрастает;

dy = = 0 при r = rэ, max y(r);

dr 0 при r rэ, y(r) убывает;

8.2. Характеристики уравнения переноса если sin c+ 0, то dy 0, y(r) убывает;

dr при cos c+ 0:

если sin c+ 0, то dy 0, y(r) возрастает;

dr если sin c+ 0, то 0 при r rэ, y(r) возрастает;

dy = = 0 при r = rэ, min y(r);

dr 0 при r rэ, y(r) убывает.

Если y(r) берется со знаком «», то при cos c+ 0:

если sin c+ 0, то dy 0, y(r) убывает;

dr если sin c+ 0, то 0 при r rэ, y(r) возрастает;

dy = = 0 при r = rэ, max y(r);

dr 0 при r rэ, y(r) убывает;

при cos c+ 0:

если sin c+ 0, то 0 при r rэ, y(r) возрастает;

dy = = 0 при r = rэ, min y(r);

dr 0 при r rэ, y(r) убывает;

+ если sin c3 0, то dy 0, y(r) возрастает.

dr Кривая y(r) достигает максимума при rэ или на ветви со знаком «+», если cos c+ 0 и sin c+ 0, или на ветви со знаком «», если cos c+ 3 3 и sin c+ 0;

в этих случаях max y(r) = c4.

r=rэ Кривая y(r) достигает минимума при rэ на ветви со знаком «+», если cos c+ 0 и sin c+ 0, или на ветви со знаком «», если cos c+ 0 и 3 3 sin c+ 0;

в этих случаях min y(r) = c4.

r=rэ Если условно называть «верхней» ветвью y(r) часть кривой, лежащей выше y = y0, а «нижней» – ниже y = y0, то для общности естественно – 574 Глава 8. Сферические задачи с осевой симметрией уравнение для «верхней» ветви брать в виде c1 c4 cos c+ + c4 | sin c+ | r2 c 3 (8.131) y(r) =, r а для «нижней»

c1 c4 cos c+ c4 | sin c+ | r2 c 3 (8.132) y(r) =.

r Оказывается, что знак sin c+ позволяет определить направление луча s. Лучи с одинаковыми cos c+ и c sin c+, отличающимися только знаками, 3 проходят одинаковые пути в области r, y, но только в противоположных направлениях.

Графики кривых при различных параметрах c+ схематически изображены на рис. 8.1 (a, b, c, d, e, f ).

Рассмотрим отдельно некоторые частные случаи.

Случай sin c+ = 0, cos c+ = 1 (cos = 0):

3 c1 c y(r) 0, если c1 = 0 или c4 = 0;

0, y(r) = r cc r(y) = 1 4 0, r(y) 0, если c1 = 0 или c4 = 0;

y sin y = r1 cos 1 ± c2 y 2 ;

(8.133) y max y = c4 достигается при r = c1 ;

асимптотой является y = 0 (рис. 8.1,g).

Случай sin c+ = 1, cos c+ = 0 (cos = c4 ):

3 ±c4 r 2 c2 c2 y c 1 y(r) =, r(y) = ;

c2 y r y как обычно (8.123) или (8.124);

y(r) экстремумов при конечных r не имеет, симметрична относительно оси r;

ее асимптоты y = ± cos = ±c4 (рис. 8.1,h).

Случай sin c+ = 0, cos c+ = 1 (cos = 0):

3 c1 c y(r) = y(r) 0, если c1 = 0 или c4 = 0;

0, r cc r(y) = 1 4 0, r(y) 0, если c1 = 0 или c4 = 0;

y y как в (11.18);

min y(r) = c4 достигается при r = c1 ;

асимптотой y(r) является y = 0 (рис. 8.1,k).

Рассмотрение случая cos 1 0 (c3 = c = const) проводится аналогично случаю cos 1 0. Исходными выражениями здесь являются (8.109), (8.111) и вместо (8.117) имеем r cos = c1 c4 cos c ± c4 sin c r 2 c2.

3 8.2. Характеристики уравнения переноса y Б sin c3 c | cos | C y D rэ r c1r | cos | sin c3 a y Б c y0 sin c3 C | cos | D c1 rэ r 0 r | cos | sin c3 b y Б c sin c3 y0 = | cos | C D r c1 r = r э | cos | sin c3 c Рис. 8. 576 Глава 8. Сферические задачи с осевой симметрией y sin c3 | cos | rэ c1 r r 0 D C y | cos | sin c3 c Б d y sin c3 | cos | rэ c1 r r D | cos | C y0 sin c3 c Б e y sin c3 | cos | c1 rэ = r r D y0 = | cos | C sin c3 c Б f Рис. 8.1. Продолжение 8.2. Характеристики уравнения переноса y Б y0 = c cos = c1 = rэ r g y c4 = | cos | D y0 = 0 r r0 = c c4 = | cos | h y c1 = rэ cos = r y0 = c Б k Рис. 8.1. Продолжение 578 Глава 8. Сферические задачи с осевой симметрией Параметрические уравнения r() и y() и выражения для r и y остаются 0, только c+ следует заменить на c.

такими же, что и в случае cos 1 3 Из уравнения (8.126) аналогично получаются уравнения луча в виде:

если y 2 = cos2, то c1 c r(y) = ;

2y cos c если y 2 = cos2, то c1 c4 cos c y ± c1 cos c2 y r(y) =, y 2 cos причем если sin c 0, «+», если sin c 0, «», иначе y cos c | sin c | c2 y 3 (8.134) r(y) = c1 c4, y 2 cos или в виде c1 c4 cos c ± cos r 2 c (8.135) y(r) = ;

r область значений r, y вдоль луча:

y2 c2.

0 c1 r +, При r = c1 cos c 0, 0, если cos c = 0, y(c1 ) = y0 = c4 cos c3 = 0, если cos c 0.

0, если При r имеются асимптоты c1 c4 cos c c ± cos 1 = ± cos ;

lim y(r) = lim r r r r c y = 0 при r 0 =.

| sin c | В случае sin c = 0 y(r) нуль не пересекает. Производная c1 c4 cos c r2 c2 ± c2 cos dy =.

dr r 2 r 2 c dy Условие = 0 – условие экстремума.

– dr Если y(r) берется со знаком «+», то должно быть c1 c4 cos c rэ c2 + c2 c4 sin c = 3 8.2. Характеристики уравнения переноса или cos c rэ c2 = c1 sin c 3 – это возможно, если sin c и cos c разных знаков.

– 3 В случае знака «» должно быть c1 c4 cos c rэ c2 c2 c4 sin c = 3 или cos c rэ c2 = c1 sin c 3 – это возможно, если sin c и cos c одинаковых знаков.

– 3 Для обоих случаев c2 c 2 или rэ = rэ =.

| cos c | cos2 c Если cos c = 0, то экстремумы достигаются только при r.

Если y(r) берется со знаком «+», то при cos c 0:

если sin c 0, то dy 0, y(r) возрастает;

dr если sin c 0, то 0 при r rэ, y(r) возрастает;

dy = = 0 при r = rэ, min y(r);

dr 0 при r rэ, y(r) убывает;

при cos c 0:

если sin c 0, то 0 при r rэ, y(r) возрастает;

dy = = 0 при r = rэ, max y(r);

dr 0 при r rэ, y(r) убывает;

если sin c 0, то dy 0, y(r) убывает.

dr Если y(r) берется со знаком «», то при cos c 0:

если sin c 0, то 0 при r rэ, y(r) возрастает;

dy = = 0 при r = rэ, min y(r);

dr 0 при r rэ, y(r) убывает;

580 Глава 8. Сферические задачи с осевой симметрией если sin c 0, то dy 0, y(r) возрастает;

dr при cos c 0:

если sin c 0, то 3 dy 0, y(r) убывает;

dr если sin c 0, то 0 при r rэ, y(r) возрастает;

dy = = 0 при r = rэ, max y(r);

dr 0 при r rэ, y(r) убывает.

Кривая y(r) достигает максимума при r = rэ или на ветви со знаком «+», если cos c 0 и sin c 0, или на ветви со знаком «», если cos c 3 3 и sin c 0;

при этом max y(r) = c4.

r=rэ Кривая y(r) достигает минимума при r = rэ или на ветви со знаком «+», если cos c 0 и sin c 0, или на ветви со знаком «», если cos c 3 3 и sin c 0;

при этом min y(r) = c4.

r=rэ Уравнение для «верхней» ветви c1 c4 cos c + c4 | sin c | r2 c 3 (8.136) y(r) =, r а для «нижней» ветви c1 c4 cos c c4 | sin c | r2 c 3 (8.137) y(r) =.

r Два луча с одинаковыми по величине и знаку cos c и с одинаковыми по величине, но отличными по знаку sin c проходят одинаковые пути в области r, y, но в противоположных направлениях.

Графики кривых те же, что и в случае cos 1 0 (рис. 8.1, a, b, c, d, e, f ).

Частные случаи.

Случай sin c = 1, cos c = 0 (cos = c4 ):

3 ±c4 r 2 c2 c2 y c 1 y(r) = ;

r(y) = ;

c2 y r y как обычно (8.123) или (8.124);

y(r) экстремумов при конечных r не имеет, симметрична относительно оси r;

ее асимптоты y = ± cos (рис. 8.1,h).

Случай sin c = 0, cos c = 1 (cos = 0):

3 c1 c y(r) 0, если c1 = 0 или c4 = 0;

0, y(r) = r cc r(y) 0, если c1 = 0 или c4 = 0;

r(y) = 1 4 0, y 8.3. Алгоритм вычисления криволинейных координат y как в (8.133);

max y = c4 достигается при r = c1 ;

«верхняя» и «нижняя»

ветвь y(r) совпадают;

ее асимптота y = 0 (рис. 8.1,g).

Случай sin c = 0, cos c = 1 (cos = 0):

3 c1 c y(r) = y(r) 0, если c1 = 0 или c4 = 0;

0, r cc r(y) = 1 4 0, r(y) 0, если c1 = 0 или c4 = 0;

y y как в (8.133);

min y(r) = c4 при r = c1 ;

асимптота y = 0 (рис. 8.1,k).

Для наглядности дадим геометрическую интерпретацию кривым (a k) рис. 8.1.

К классу лучей, имеющих траекторию вида рис. 8.1,a, b, c, относятся лучи, у которых точка C наименьшей удаленности от центра O (ей соответствует r = c1 ) лежит в верхней полусфере /2, а у лучей с траекториями рис. 8.1,d, e, f точка C находится в нижней полусфере /2.

Лучи с траекториями вида рис. 8.1,h располагаются в плоскостях, парал лельных оси z;

точка C в этом случае находится на линии = /2.

Траектории вида рис. 8.1,g, k соответствуют лучам, проходящим в плос костях, перпендикулярных оси z, но пересекающих ось z выше или ниже центра O. При этом лучи рис. 8.1,g находятся в верхней полусфере /2, а лучи рис. 8.1,k – в нижней полусфере /2.

– Особый класс составляют лучи, траектория которых в плоскостях r, y вырождается в прямую y = 0. Им соответствует cos = 0. Такие лучи лежат в плоскости Qhor, проходящей через центр O перпендикулярно оси z.

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

§ 8.3. Алгоритм вычисления криволинейных координат на траекториях характеристик Рассматривается задача переноса солнечного излучения в земной атмосфере в масштабах всей планеты. В настоящем разделе ограничиваемся приближением сферической модели с осевой симметрией, когда пространственные координаты и направления распространения излучения описываются сферическими систе мами координат. При смещении вдоль лучей, совпадающих с характеристиками дифференциального оператора переноса, изменяются, как правило, все четыре 582 Глава 8. Сферические задачи с осевой симметрией переменные. Далее описан алгоритм расчета криволинейных координат на траекториях характеристик в четырех-мерном фазовом пространстве.

8.3.1. Алгоритм вычисления, r, y,, t. Прежде всего введем некоторые обозначения для разностной сети задачи:

l = 0 l0 L, L = 1;

0 = 1, l0 = 0, l = cos l, i = I 0 I, I = 1, 0 = 0, I = 1;

i = cos i, j = J 0 J, tJ = 1, t0 = 0, tJ = 1.

tj = cos j, Далее будем считать, что координатами рассматриваемой точки А являются y1 = l и r1 = rk, а направление рассматриваемого луча s, проходящего через точку A(r1, y1 ), определяется cos 1 = i и cos 1 = tj. По этим данным с помощью результатов, полученных выше, можно получить первые интегралы оператора переноса и следующую информацию для луча s:

1 cos2 1 = r1 sin 1 ;

c1 = r 1 cos2 1 1 y1 ;

c2 = sin 1 sin 1 = c2 y y1 1 c2 ;

c4 = cos x = ;

sin x = ;

2 c4 c y 0 x = arccos ;

c если cos 1 0, то находим c+ = x 1, 3 cos c+ = sin (x + 1 ) = sin x cos 1 + cos x sin 1, sin c+ = cos (x + 1 ) = cos x cos 1 sin x sin 1 ;

8.3. Алгоритм вычисления криволинейных координат если cos 1 0, то находим c = x + 1, 3 cos c3 = sin (x 1 ) = sin x cos 1 cos x sin 1, sin c = cos (x 1 ) = cos x cos 1 + sin x sin 1 ;

c4 sin c+, если cos 1 0, cos = c4 sin c, если cos 1 0;

c4 cos c+, если cos 1 0, y0 = c4 cos c, если cos 1 0;

c, если cos 1 0, | cos c+ | rэ = c1, если cos 0;

| cos c | c, если cos 1 0, | sin c+ | r0 = c1, если cos 0.

| sin c3 | При анализе прохождения луча через ячейку будем считать, что эта информация нам уже известна. Через r, y,, t мы обозначаем значения аргументов функции gr на границе, через которую характеристика входит в ячейку, причем y = cos, = cos, t = cos. В удобном для нас виде выше выписаны основные формулы, которые будут использованы в алгоритме вычисления, r, y,, t под нумерацией, установленной ниже.

8.3.2. Формулы для вычисления y. В дальнейшем будут использоваться следующие формулы для вычисления y :

если |y| = | cos |, то (формула (y 1)) y 2 cos 1 y1 cos ± |y sin 1 | c2 y (8.138) y = r1, y 2 cos где выбор знака определяется признаком 0, если « + », = если « »;

1, если |y| = | cos |, то (формула (y 2)) r1 (cos2 y1 ) r1 (cos2 y1 ) 2 (8.139) y = = ;

2 cos (cos 1 cos y1 ) 2c4 cos sin 1 cos c (y 3) – упрощенная форма для (y 2) (8.139), если y1 = 0:

– r (8.140) y = ;

2 cos 584 Глава 8. Сферические задачи с осевой симметрией (y 4) – упрощенная форма для (y 1) (8.138), если y = 0:

– r1 cos (8.141) y = ;

cos (y 5) – эквивалентная форма для (y 1) (8.138), если cos = 0:

– sin y = r1 cos 1 ± c2 y 2 (8.142), y 0, если « + », = если « »;

1, в таких ситуациях, как правило, y = 0.

8.3.3. О выборе. Если |y| = | cos |, то всегда имеются два y, соответствующих одному значению |y|, так как y зависит только от |y|;

причем имеются две возможности: либо y2 y1 0, либо y1 0, y2 0.

Мы всегда будем выбирать y = y1 0.

При |y| | cos | два значения y получаются при одном и том же y;

при |y| | cos | одно y соответствует y, а другое соответствует y (см. графики a–k рис. 8.1).

Если |y| = | cos |, то возможны три значения y, два из которых равны ± и лишь одно конечно, причем последнее больше нуля;

оно и принимается за y (см. формулу (y 2)).

Анализируя формулу (y 1) легко установить, что если |y| | cos |, то если 0 y1 y2, то y1 соответствует = 1, если y1 0, y2 0, то y1 соответствует = 0;

если |y| | cos |, то если 0 y1 y2, то y1 соответствует = 0, если y1 0, y2 0, то y1 соответствует = 1.

8.3.4. Формулы для вычисления r. Будем пользоваться в дальнейшем следующими формулами для вычисления r :

если cos 1 = 1, то (формула (r 1)) r = r r1 ;

(8.143) если 1 cos 1 1, то (формула (r 2)) r = r1 cos 1 ± r 2 c2, (8.144) а знак выбирается с помощью признака 0, если « + », r = если « »;

1, 8.3. Алгоритм вычисления криволинейных координат если cos 1 = 1, то (формула (r 3)) r = r1 r. (8.145) 8.3.5. Формулы для вычисления r(), y(), (r), t(y). Обозначим через (1 ) формулу r() = 2 2b + c, 2, (8.146) b = r1 cos 1, c = r r [rk, rk +1 ];

через (2 ) – формулу – r y cos y() = 1 1, (8.147) r() y [, 1 ];

через (3 ) – формулу – r 2 c, (r) = ± r (8.148) 0, если « + », 1, если « », = [, +1 ];

через (4 ) – формулу – c2 y 2 t(y) = ± y = 1,, 1 y2 (8.149) 0, если « + », t = если « », 1, t [tx, tx+1 ].

8.3.6. О выборе r, и r. Возможны следующие ситуации:

(1) если 0 cos 1 1 и c1 rk1, то, так как (r) 0, полагается r = 1, = 0;



Pages:     | 1 |   ...   | 10 | 11 || 13 | 14 |   ...   | 15 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.