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

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

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


Pages:     | 1 |   ...   | 4 | 5 || 7 | 8 |   ...   | 15 |

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

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

§ 4.3. Интегрирование по характеристике комплексного уравнения переноса Краевые задачи для ПЧХ формально совпадают с одномерной краевой задачей для уравнения переноса излучения в плоском слое с тем лишь принципиальным различием, что искомая функция (z, p, s) и коэффици ент экстинкции (p, z) t (z) i(p, s ) являются комплекснозначными и зависящими от действительного параметра p = {px, py }. Для численного решения краевых задач для ПЧХ при фиксированных значениях парамет ров – пространственных частот – используем принцип построения алгоритма – – решения скалярного уравнения переноса, известный как метод итераций с интегрированием уравнения по характеристике дифференциального оператора переноса и квадратурами на единичной сфере.

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

Комплексные уравнения переноса для ПЧХ имеют одинаковую структуру:

выделен одномерный интегродифференциальный оператор L(p) Lz (p) s (z)M, коэффициенты которого не зависят от горизонтальной координаты r, а от пространственной частоты p зависят как от параметра. Это предопределяет общие черты алгоритма их решения.

В каждой из ПЧХ – решении задачи вида – L(p)n (z, p, p, s) = Fn1 (z, p, p, s), (4.16) n = 0, n = q RH n + fn1 (p, p, s), 0 H где через Fn1, fn1 обозначены источники, описывающие рекуррентную связь между ПЧХ двух последовательных приближений, а p – набор параметров, – 4.3. Интегрирование по характеристике комплексного уравнения отличных от параметра p, входящего в оператор переноса Lz (p), выделяются две составляющие:

n = 0 + g.

n n Интегрирование задачи для нерассеянной компоненты Lz (p)0 = Fn1 (z, p, p, s), 0 = 0, (4.17) = fn1 (p, p, s) n n n 0 H осуществляется по характеристикам оператора Lz (p) c учетом конкретной структуры источников Fn1, fn1 по переменным z,, и зависимости от параметров p, p. Для решения задач обычно разрабатывается отдельный мо дуль в программах, вычисляющий решение либо по квадратурным формулам, либо по аналитическим выражениям, которые предварительно выводятся.

Задачи для определения многократно рассеянной, диффузной составляю щей ПЧХ L(p)g = S0, = q RH g + RH g g = 0, (4.18) n n n n n n 0 H оказываются одинаковыми при различных источниках в (4.16), если интегралы столкновений S0 и граничные значения RH 0 задаются в дискретном, n n сеточном виде.

Вводится разностная сеть задачи = {zk, j, i } и определяются сеточные функции zk b± sin ± (px cos i + py k 1 = 0, sin i ), t (z) dz, K = (H), ji j 0± 0 (zk, ±, i ), g± g (zk, ±, i ).

n n j ji kji kji Коэффициенты задачи t (z), s (z), a (z), g(z), (z, ) считаются кусочно постоянными по высоте, т. е. для z [zk, zk+1 ] полагается g(z) = gk, t (z) = tk, s (z) = sk, a (z) = ak, для z [zl, zl+1 ] (z, ) = l ().

Обращение дифференциального оператора в уравнениях (4.17), (4.18) осу ществляется путем интегрирования по характеристикам, которые совпадают с траекториями лучей s = s± = {±, i }. Для сокращения записей индексы ij j луча i, j в сеточных функциях опускаем.

Обозначим B Sg + S0, BH q Rg + R n n n n и представим каждую из комплекснозначных функций = g, B, Bn в виде n действительной (Re) и мнимой (Im) компонент, принимающих действительные значения:

= Re + iIm, B = B Re + iB Im. (4.19) Интегрированием уравнения (4.18) по характеристике дифференциального оператора Lz (p), определяемой условиями = const, = const, получаем 268 Глава 4. Комплексное уравнение переноса основные формулы метода характеристик:

z + B(t, +, ) exp ((t) (z)) /+ dt, (z,, ) = + H (z,, ) = B(t,, ) exp ((z) (t)) /| | dt + (4.20) | | z + BH (, ) exp ((z) (H)) /| |, где комплексная «оптическая толщина»

z (z) (t) dt = (z) ibz.

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

+ = 0, k=1: 1ji zk+ + (z ) dz + k = 1 (K 1) : + + k+1,ji = kji exp + j zk (4.21) zk+1 zk+ 1 Bji (t) exp + + (t ) dt dt, + + + j j zk t = BHji, k=K : Kji zk+ (z ) dz + k = (K 1) 1 : = k+1,ji exp | | kji j (4.22) zk zk+1 t 1 Bji (t) exp (t ) dt dt.

+ | | |j | j zk zk В этих формулах введено ± (z) t (z) ib±.

ji Учитывая кусочно-постоянную аппроксимацию по высоте коэффициента экстинкции t (z), выделим в экспоненциальных множителях действительные 4.3. Интегрирование по характеристике комплексного уравнения и мнимые части:

zk+ exp + + (z ) dz = exp k /+ cos d+ + i sin d+, (4.23) j kji kji j zk zk+ exp + (t ) dt = exp tk (t zk+1 )/+ + + cos gk+1,ji + i sin gk+1,ji, + j j t d+ b+ zk /+ = ji zk, + ji b+ /+ = tg + (px cos i + py sin i ), + ji j ji j j kji gk+1,ji ji (zk+1 t) = b+ (zk+1 t)/+.

+ + ji j Подставим комплексные выражения для функций В, (4.19) и экспонен циальных множителей (4.23) в формулы (4.21), (4.22) и, приравняв отдельно мнимые и действительные члены в обеих частях равенства, получим:

Re+ = Re+ cosd+ Im+ sind+ exp k /+ + j k+1,ji kji kji kji kji zk+ Bji (t)cosgk+1,ji Bji (t)singk+1,ji exp tk (tzk+1 )/+ dt, + + Re+ Im+ ++ j j zk Im+ = Re+ sind+ +Im+ cosd+ exp k /+ + j k+1,ji kji kji kji kji zk+ Bji (t)singk+1,ji +Bji (t)cosgk+1,ji exp tk (tzk+1 )/+ dt. (4.24) + + Re+ Im+ ++ j j zk Аналогичные преобразования приводят к выражениям для направлений s :

Re = Re cos d Im sin d exp k /| | + j k,ji k+1,ji kji k+1,ji kji zk+ Bji (t) cos gkji Bji (t) sin gkji exp tk (zk t) /| | dt, Re Im + j |j | zk Im = Re sin d + Im cos d exp k /| | + j kji k+1,ji kji k+1,ji kji zk+ Bji (t) sin gkji + Bji (t) cos gkji exp tk (zk t) /| | dt, Re Im + j |j | (4.25) zk 270 Глава 4. Комплексное уравнение переноса где введены следующие представления:

zk+ = exp k /| | cos d + i sin d, exp (z ) dz j kji kji |j | zk t = exp tk (zk t)/| | exp (z ) dz cos gkji + i sin gkji, | | j j zk dkji bji zk /| | = ji zk ji b /| | = | tg | (px cos i + py sin i ),, j ji j j gkji ji (t zk ) = b (t zk )/| |.

(4.26) ji j Конечно-разностный аналог для формул метода характеристик (4.24), (4.25) зависит от способа аппроксимации функции источника B(z) на [zk, zk+1 ]. Это может быть кусочно-постоянная, кусочно-линейная, кусочно экспоненциальная, дробно-рациональная или еще какая-то аппроксимация.

В случае кусочно-постоянной аппроксимации высотной зависимости функ ции источника на [zk, zk+1 ]:

Bji ± (z) = Bk+1,ji + Bkji± /2, Re ± Bji ± (z) = Bk+1,ji + Bkji± / Im ± Re Re Im Im получаем следующую двухточечную схему:

Re + = 0, Im + = 0 ;

k = 1: 1ji 1ji k = 1 (K 1): Re + = Re + Dkji Im + Ekji + + + k+1,ji kji kji + Bk+1,ji + Bkji+ A+ Bk+1,ji + Bkji+ Ckji, Re + Im + + Re Im kji Im + = Im + Ekji + Im + Dkji + + + k+1,ji kji kji + Bk+1,ji + Bkji+ Ckji + Bk+1,ji + Bkji+ A+ ;

Re + + Im + Re Im kji Re = BHji, Im = BHji ;

Re Im k = K: Kji Kji k = (K 1) 1: Re = Re Dkji Im Ekji + kji k+1,ji k+1,ji + Bk+1,ji + Bkji A Bk+1,ji + Bkji Ckji, Re Im Re Im kji Im = Im Ekji + Im Dkji + kji k+1,ji k+1,ji + Bk+1,ji + Bkji Ckji + Bk+1,ji + Bkji A. (4.27) Re Im Re Im kji Коэффициенты A±, Dkji, Ckji, Ekji представляются аналитически и ± ± ± kji зависят от коэффициентов уравнения переноса, параметра p и разностной 4.3. Интегрирование по характеристике комплексного уравнения сети задачи:

Dkji = cos d+ exp k /+, + Ekji = sin d+ exp k /+ +, j j kji kji zk+ A+ = cos gk+1,ji exp tk (t zk+1 )/+ dt = + 2+ j kji j zk tk cos d+ b+ sin d+ exp k /+ / 2 tk + (b+ ) = tk, ji j ji kji kji zk+ + sin gk+1,ji exp tk (t zk+1 )/+ dt = + Ckji = 2+ j j zk b+ tk sin d+ + bji cos d+ exp k /+ / 2 tk + (b+ ) =, ji j ji kji kji Dkji = cos d exp k /| |, Ekji = sin d exp k /| |, j j kji kji zk+ A = cos gkji exp tk (zk t)/| | dt = 2| | j kji j zk tk cos d b sin d exp k /| | / 2 tk + (b ) = tk, ji j ji kji kji zk+ sin gkji exp tk (zk t)/| | dt = Ckji = 2| | j j zk = b tk sin d + b cos d exp k /| | /[2(tk + (b )2 )]. (4.28) ji ji j ji kji kji Если на [zk, zk+1 ] ввести линейную аппроксимацию высотной зависимости функции источника:

Bji ± (z) = Bkji± + Bk+1,ji Bkji± (z zk )/zk, Re ± Re Re Re Bji ± (z) = Bkji± + Bk+1,ji Bkji± (z zk )/zk, Im ± Im Im Im то конечно-разностный аналог формулы интегрирования по характеристике приводится к следующему виду:

Re + = Re + Dkji Im + Ekji + Bk+1,j + Bkji+ A+ Bk+1,ji + Bkji+ Ckji + + + Re + Im + + Re Im k+1,ji kji kji kji A+ zk A+ + + Ckji zk Ckji Bk+1,ji Bkji+ Re + Bk+1,j Bkji+ ;

Im + kji kji Re Im + zk zk Im + = Re + Ekji + Im + Dkji + Bk+1,ji + Bkji+ Ckji + Bk+1,ji + Bkji+ A+ + + + Re + + Im + Re Im k+1,ji kji kji kji + + A+ zk A+ Ckji zk Ckji Bk+1,ji Bkji+ + Re + Bk+1,ji Bkji+, (4.29) Im + kji kji Re Im + zk zk 272 Глава 4. Комплексное уравнение переноса где коэффициенты D +, E +, A+, C + совпадают с (4.28), a для коэффициентов A+, C + получаем аналитические выражения zk+ A+ (zk+1 t) cos ji (zk+1 t) exp tk (t zk+1 )/+ dt = + =+ j kji zk = + / tk + (b+ )2 + k /+ + + cos d+ j ji j kji kji kji (d+ + kji ) sin d+ exp k /+ + (4.30), j kji kji zk+ + (zk+1 t) sin ji (zk+1 t) exp tk (t zk+1 )/+ dt = + Ckji = + j j zk = + / tk + (b+ )2 + d+ + kji cos d+ + + kji j ji kji kji + k /+ + + sin d+ exp k /+, j j kji kji где + tk (b+ )2 / tk + (b+ )2, kji 2tk b+ / tk + (b+ )2.

+ 2 2 ji ji ji ji kji Для направлений s :

Re =Re Dkji Im Ekji + Bk+1,ji +Bkji A Bk+1,ji +Bkji Ckji Re Re Im Im kji k+1,ji k+1,ji kji A zk A Ckji zk Ckji kji kji Bk+1,ji Bkji + Bk+1,ji Bkji, Re Re Im Im zk zk Im =Re Ekji +Im Dkji + Bk+1,ji +Bkji Ckji + Bk+1,ji +Bkji A Re Re Im Im kji k+1,ji k+1,ji kji A zk A Ckji zk Ckji kji kji Bk+1,ji Bkji Bk+1,ji Bkji.

Re Re Im Im (4.31) zk zk Коэффициенты A, E, D, C совпадают с (4.28), а для коэффициентов C, A получаем zk+ A = cos ji (t zk ) exp tk (zk t)/| | dt = j kji |j | zk | | k /| | + cos d = j j kji kji kji d + kji sin d exp k /| | / tk + (b )2, 2 (4.32) j ji kji kji 4.3. Интегрирование по характеристике комплексного уравнения zk+ (t zk ) sin ji (t zk ) exp tk (zk t)/| | dt = Ckji = j |j | zk = | | kji d + kji cos d + j kji kji + k /| | + sin d exp k /| | / tk + (b )2, j j ji kji kji где tk (b )2 / tk + (b )2, kji 2tk b / tk + (b )2.

+ 2 2 ji ji ji ji kji В случае однородного слоя по высоте задачу (4.16) можно записать для пространственной координаты – оптической толщины по вертикали, положив – p = {px /t, py /t }, z =, t = 1, s – альбедо акта рассеяния. Если ввести – равномерный шаг по высоте zk = const и обозначить cos+ cos ji, + ji sin+ sin ji, то формулы (4.28) для коэффициентов разностной схемы + ji существенно упрощаются:

Dji = cos+ exp /+, + Eji = sin+ exp /+, + ji j ji j A+ = 1 cos+ b+ sin+ exp /+ / 2 1 + (b+ )2, ji ji ji ji j ji Cji = b+ sin+ +b+ cos+ exp /+ + / 2 1 + (b+ )2.

ji ji ji ji j ji Упростятся и формулы (4.30):

A+ = + + /+ + + cos+ ji j ji j ji ji (ji + ji ) sin+ exp /+ + + 1 + (b+ )2, ji j ji Cji = + ji + + ji + ji cos+ + + + j ji + /+ + + sin+ exp /+ / 1 + (b+ )2, j ji ji j ji где + 1 (b+ )2 / 1 + (b+ )2, ji 2b+ / 1 + (b+ )2.

+ ji ji ji ji ji Аналогичны упрощения для s :

Dji = cos exp /| |, Eji = sin exp /| |, ji j ji j A cos b sin /| | (b )2 ) = 1 / 2(1 + exp, ji ji ji ji j ji b sin b cos /| | (b ) / 2 1+ Cji = exp, ji ji ji ji j ji A = | | /| | + cos ji j ji j ji ji ji + ji sin exp /| | / 1 + (b )2, ji j ji Cji = | | ji ji + ji cos + j ji + /| | + sin exp /| | / 1 + (b )2, j ji ji j ji 274 Глава 4. Комплексное уравнение переноса где 1 (b )2 2b + ji ji ji,.

1 + (b )2 1 + (b ) ji ji ji Если ввести дополнительное ограничение, взяв симметричную сеть по углу такую, что + = | |, j = 1 J, J = J + = J (такая ситуация обычна j j при расчетах с гауссовой сетью j одинакового порядка для интервалов [0, 1] и [1, 0]), и обозначить sin j = sin + = sin, j = + = | |, j j j j cosji = cos+ = cos, sinji = sin+ = sin, ji ji ji ji bji = b+ = b, + ji = ji = ji, j = exp (/j ), ji ji то получим + + Dji = Dkji = Dkji = j cosji, Eji = Ekji = Ekji = j sinji, Aji = A+ = A = [1 j (cosji bji sinji )] / 2 1 + b2, ji kji kji + Cji = Ckji = Ckji = [bji j (sinji +bji cosji )] 2 1 + b2, ji j Aji = A+ = A = {ji j [(/j + ji ) cosji (ji + ji ) sinji ]}, kji kji 1 + b ji j + {ji j [(ji + ji ) cosji + (/ji + ji ) sinji ]}, Cji = Ckji = Ckji = 1 + b ji где 1 b2 2bji ji = + = = + ji, ji = ji = ji =.

ji ji b2 1 + b 1+ ji ji Если px = 0, py = 0, то приходим к обычной плоской задаче с коэффициентами bji = 0, cosji = 1, sinji = 0, Dji = j, Aji = (1 j ) /2, Eji = 0, Cji = и конечно-разностному аналогу:

Re + = Re + j + Bk+1,ji + Bkji+ Aji, Re + Re k+1,ji kji Re = Re j + Bk+1,ji + Bkji Aji, Re Im + = Im = 0.

Re kji k+1,ji kji kji Расчетные формулы полностью совпадают с решением обычного уравнения переноса для плоского слоя (см. гл. 1). Обращаем внимание на тот факт, что в конечно-разностной схеме интегрирования комплексного уравнения переноса по характеристике переход к значениям параметров px = 0, py = 0 не приводит к авостным ситуациям, поэтому в этом частном случае можно использовать общие расчетные формулы.

4.3. Интегрирование по характеристике комплексного уравнения Решение задачи (4.18) проводится методом итераций – последовательных – (N ) – приближенное решение, приближений по кратности рассеяния. Если – полученное на итерации с номером N, то следующее приближение находится по схеме (N ) Lz (N +1) = B (N ), (N +1) (N +1) = 0, (4.33) = BH, 0 H в которой функция источника (N ) B (N ) S(N ) + S0, = q R(N ) + R BH.

Отметим важное достоинство предложенного алгоритма. Построенные конечно-разностные аналоги метода характеристик сохраняют последователь ность приближений и не перемешивают итерации: все расчетные формулы имеют вид рекуррентных соотношений, когда значения сеточных функций Re ±, Im ± на следующем k-м уровне слоя на итерации с номером N kji kji вычисляются через значения функций Re ±, Im ±, взятых с предыдущего k±1ji k±1ji (k ± 1)-го уровня слоя на той же итерации. Произошло четкое разделение действительных и мнимых компонент решения, не нарушающее монотонности характеристической схемы.

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

Но компоненты Re, Im оказываются колеблющимися по, функциями, и чем больше значения px, py, тем больше колебаний. Поэтому в качестве критерия сходимости итераций проверяется условие (N +1)± (N )± (N +1)± ± = |kji | |kji | /|kji |, max {zk,±,i } j где – заданная относительная точность итераций, для слабо колеблющейся – функции 2 1/ ± = Re ± + Im ±.

kji kji kji С увеличением оптической толщины слоя и уменьшением величины px, py сходимость итераций ухудшается. В этих случаях рекомендуется воспользоваться ускорением сходимости итераций в виде верхней релаксации, когда последовательные приближения вводятся по схеме (4.33), а интегралы (N ) B (N ), BH вычисляются через функции (N ) = (N ) + (N ) (N 1), где релаксационный параметр принимает значения в диапазоне 1 2.

Из анализа формул для коэффициентов разностного аналога устанавли вается, что ± bji A±, ± Ckji, + (b± ) kji tk ji 276 Глава 4. Комплексное уравнение переноса т. е. lim A± = 0, lim Ckji = 0 при p ±, и при достаточно больших ± kji значениях px, py вклад кратного рассеяния, который описывается функциями B Re ±, B Im ±, становится малым.

В комплексных уравнениях переноса для ПЧХ пространственная пере менная z является геометрической величиной, измеряемой в единицах длины (для задач атмосферной оптики обычно это километры). В задачах теории переноса для плоского слоя часто от геометрической переменной переходят к оптической толщине (z). Сделаем такой переход в задаче (4.15), поделив все члены уравнения на величину t (z) и использовав соотношение d = t (z) dz.

В результате получим уравнение (4.34) + (, p)(, p, s) = ( ) M + c( )V1 (p) + + ( ) M 00 00 + c( ), где альбедо акта рассеяния s ( (z)) g( (z)) ( ) =, ( ) =, t ( (z)) t ( (z)) а коэффициент экстинкции (p, s ) b (, p) = 1 i =1i.

t ( (z)) t ( (z)) В однородном по высоте слое t (z) = t = const, s (z) = s = const и, следовательно, ( ) = = const, ( ) = = const, (, p) = (p).

Приняв p = {x, py }, где px = px /t, py = py /t (параметры px, py p безразмерные, так как исходные параметры px, py имеют ту же размерность, что и коэффициент экстинкции t (z)), можно использовать без изменений алгоритм, описанный выше, если в расчетных формулах положить tk = 1, sk =, gk =, zk = k.

В случае неоднородного слоя получается алгоритм, аналогичный изло женному выше, если ввести кусочно-постоянную аппроксимацию для t ( ) на интервале [k, k+1 ]: t ( ) = tk.

Основные формулы интегрирования комплексного уравнения переноса по характеристике примут следующий вид:

1 + + + (t ) dt B(t,, ) exp + (,, ) = + dt t – аналог формулы (4.20) для s +, – (H) t 1 (,, ) = (t ) dt B(t,, ) exp dt + | | | | (H) + BH (, ) exp (t) dt | | 4.3. Интегрирование по характеристике комплексного уравнения – аналог формулы (4.20) для s, ± ( ) 1 ib± /t ( ). Эти выражения – приводятся к дискретным аналогам для разностной сети k с шагом k = k+1 k подобно тому, как это изложено выше.

Сопоставляя конечно-разностные схемы, представленные через геометри ческую переменную z и оптическую толщину, находим, что все формулы совпадают, за исключением выражений для коэффициентов A±, Ckji : для ± kji схемы с перед этими коэффициентами в отличие от схемы с z стоит множитель tk. Появление этого множителя легко понять. Действительно, в задаче с z интегралы столкновений B = s M имеют размерность, равную размерности, умноженной на размерность коэффициента рассеяния s. В задаче с интегралы столкновений B = M, где безразмерная, имеют размерность, совпадающую с размерностью. Тем не менее размер ность в обеих задачах одинакова. Поэтому появившийся множитель tk в коэффициентах A±, Ckji, стоящих в разностной схеме при интегралах ± kji столкновения, приводит выражения типа B A или B C к нужной размерности (произведение tk k дает в результате sk, как в задаче с z).

Однако необходимо обратить внимание на следующую деталь. При решении задач с пространственной переменной z и оптической толщиной может получиться расхождение в результатах из-за неодинакового усреднения коэффициента экстинкции t (z) и t ( ) на шаге разностной сети zk или k соответственно. Хорошо известно, что в земной атмосфере t (z) сильно меняется по высоте (на несколько порядков). Обычно допустимо усреднение t (z) в пределах полукилометрового или километрового слоя по высоте, т. е. zk 0, 5 1 км. Соответствующие шаги k = tk zk будут сильно различаться, отвечая изменениям t (z). При этом желание уменьшить размерность K сети k или сделать сеть k более равномерной может привести к тому, что на один шаг k будет приходиться несколько шагов zk и усредненное на k значение tk будет сильно отличаться от среднего значения tk на шаге zk. Такое различие наиболее очевидно, например, в случае выбора равномерных шагов zk и k. Другое различие возникает из-за отличия аппроксимации интегралов B(z) и B( ) на интервалах [zk, zk+1 ] и [k, k+1 ].

И еще одно обстоятельство. В отличие от классической одномерной плоской задачи теории переноса, где в случае переменной достаточно знать и ( ), для решения задачи (4.34), приведенной к оптической толщине, кроме значений и ( ) требуется задание распределения t по.

Характер высотной неоднородности оказывает непосредственное влияние на получаемое решение. Этот вывод подтвержден численными расчетами. Если пытаться формально избавиться от этой зависимости, введя пространственную частоту p( ) = p/t ( ), то такая попытка не приводит к упрощению ситуации.

Напротив, возникает усложнение, связанное с тем, что на разных высотных уровнях придется иметь дело с разными пространственными частотами.

278 Глава 4. Комплексное уравнение переноса В итоге осложняются алгоритм и интерпретация расчетных результатов.

Исключение составляют задачи для однородного слоя.

§ 4.4. Расчет нерассеянной компоненты В качестве иллюстрации подхода кратко опишем алгоритмы решения краевых задач для линейных ПЧХ W10 (z, p, s) W10 + (p, z)W = (z)M W, s 10 z (4.35) W 10 = 0, W10 = q RW10 + fH 0 H и W01 (z, p, s) W01 + (p, z)W = (z)M W + (z)c(z)V (z, p) + s s 01 z + g(z) M E 00 + c(z), (4.36) W01 = 0, W01 = q RW10 + dV1 (H, p), 0 H где z s +, g(t) exp [i (p, s0 ) (z t)/0 ] dt, V1 (z, p, s0 ) = s.

0, Численное решение задачи (4.36) проводится в три этапа.

1. Находим 00 (z, s) из краевой задачи (3.39). Для нелинейных приближе ний источник в уравнении (3.40) вычисляется через предыдущее приближение.

2. При известном значении 00 определяем нерассеянную компоненту из краевой задачи Lz 0 = s (z)c(z)V1 (z, p) + g(z) M 00 00 + c(z), (4.37) 0 = 0, = q dV1 (H, p).

0 H 3. Подставляем 0 в источники задачи (4.36) и итерационным методом характеристик рассчитываем диффузную компоненту g.

При фиксированных направлениях s +, s интегрирование по характеристике уравнения (4.37) и переход к разностной сетке zk приводят 4.4. Расчет нерассеянной компоненты к следующим выражениям:

0+ = 0;

0 = f ;

k = 1: k=K:

1 K ib+ zk k /+ + k = 1 (K 1): 0+ = 0+ exp k+1 k zk+ P + (t) exp tk ib+ (t zk+1 )/+ dt;

++ (4.38) zk ib zk k /| | + k = (K 1) 1: 0 = 0 exp k k+ zk+ P (t) exp tk ib (zk t)/| | dt, + | | (4.39) zk где zk = zk+1 zk, k = k+1 k = tk zk, P ± (t) – значение правой – части, f – значение граничного условия в (4.37).

– В комплекснозначных функциях выделяем действительные и мнимые компоненты:

0± = kji + i0Im±.

0Re± kji kji В результате приходим к явным двухточечным схемам. Для фиксированного направления sij + :

0Re+ 0Im+ k = 1: = 0, = 0;

1 0Re+ k+1 = cos (b+ zk /+ ) 0Im+ sin (b+ zk /+ ) 0Re+ = k k exp (k /+ ) + Ik, Re+ k = 2 K:

0Im+ k+1 = sin (b+ zk /+ ) + k cos (b+ zk /+ ) 0Re+ 0Im+ = k exp (k /+ ) + Ik Im+. (4.40) 280 Глава 4. Комплексное уравнение переноса Для фиксированного направления sij :

0Re 0Im = f Re, = f Im ;

k=K : K K k = (K 1) 1 : k 0Re = = k+1 cos (b zk /| |) 0Im sin (b zk /| |) 0Re k+ exp (k /| |) + Ik, Re 0Im k = = k+1 sin (b zk /| |) + k+1 cos (b zk /| |) 0Re 0Im exp (k /| |) + Ik Im. (4.41) Re± Im± Через Ik и Ik обозначены соответствующие значения интегральных выражений, входящих в формулы (4.38), (4.39). Рекуррентные соотноше ния (4.40), (4.41) представляют собой конечно-разностный аналог метода характеристик для расчета нерассеянной компоненты.

Для направлений s + источник P + (t) = s (t)c(t)V1 (t) + g(t) Q+ (t) + c(t), где Q+ (t) = M E +, (t) = tk (t zk ) + k. Введем обозначение b0 (p, s0 ) = sin 0 (px cos 0 + py sin 0 ). Если 0 = 0, то b0 = px sin 0.

Рассмотрим отдельно случаи с b0 = 0 и b0 = 0.

Из формулы для V1 (t), t [zk, zk+1 ], получаем V1 (z1 ) = 0, V1 (t) = V1 (zk ) exp [ib0 (t zk )/0 ] (igk /b0 ) {1 exp [ib0 (t zk )/0 ]}. (4.42) Далее сокращаем запись: Vk V1 (zk ). Перепишем (4.38), подставляя выра жение для P + (t) и учитывая кусочно-постоянную зависимость коэффициентов 4.4. Расчет нерассеянной компоненты t (z), s (z), g(z):

k = 1 (K 1): 0+ = 0+ exp ib+ zk k /+ + k+1 k zk+ {sk S k (0 )exp[(k + sk (t zk ))/0 ] ++ zk {Vk exp[ib0 (t zk )/0 ] (igk /b0 )(1 exp[ib0 (t zk )/0 ])} dt + + gk S k (0 )exp[( + sk (t zk ))/0 ] exp tk ib+ (zk+1 t)/+ dt + zk+ gk Q+ (t)exp tk ib+ (zk+1 t)/+ dt.

++ (4.43) zk Интеграл с функцией Q+ (t) вычисляем по квадратурным формулам, а остальные интегралы можно упростить с учетом кусочно-постоянной зависи мости подынтегральных функций. Выпишем каждое слагаемое в интеграле (4.43) отдельно:

+ + + + + I + = I1 + I2 + I3 + I4 + I5.

Первое слагаемое при 0 = + и b0 = b+ I1 = C1 exp zk (sk ib0 ) /0 zk+1 sk ib+ /+ Vk + + zk+ tk ib+ /+ (tk ib0 ) / Re+ Im+ exp t dt = I1 + iI1, zk + C1 sk S k (0 ) exp (k /0 ) /+, + C1 VkRe R11 R21 VkRe R11 R Re+ Re Re Im Im I VkIm R11 R21 VkIm R11 R21, Re Im Im Re (4.44) + C1 VkRe R11 R21 + VkRe R11 R21 + Im+ Re Im Im Re I + VkIm R11 R21 VkIm R11 R21.

Re Re Im Im 282 Глава 4. Комплексное уравнение переноса Из формулы (4.42) получаем V1Re = V1Im = 0, для k = 2 K VkRe = Vk1 cos (b0 zk1 /0 ) Vk1 sin (b0 zk1 /0 ) Re Im (gk1 /b0 ) sin (b0 zk1 /0 ), (4.45) VkIm = Vk1 cos (b0 zk1 /0 ) + Vk1 sin (b0 zk1 /0 ) Im Re (gk1 /b0 ) [1 cos (b0 zk1 /0 )].

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

R11 = tk 1/+ 1/0 / tk 1/+ 1/0 + b+ /+ b0 / 2 Re, R11 = b+ /+ b0 /0 / tk 1/+ 1/0 + b+ /+ b0 / 2 Im (4.46), R21 = cos (b0 zk /0 ) exp (k /0 ) cos b+ zk /+ exp k /+, Re R21 = sin (b0 zk /0 ) exp (k /0 ) sin b+ zk /+ exp k /+.

Im Если 0 = + и b0 = b+, то R11 = + /(b+ b0 ), Re Im R11 = 0, R21 = cos (b0 zk /0 ) cos b+ zk / Re exp (k /0 ), R21 = sin (b0 zk /0 ) sin b+ zk / Im exp (k /0 ).

Если 0 = + и b0 = b+, то Re Im R11 = 0, R11 = zk, Re R21 = cos (b0 zk /0 ) exp (k /0 ), Im R21 = sin (b0 zk /0 ) exp (k /0 ).

4.4. Расчет нерассеянной компоненты Второе слагаемое при 0 = + + + I2 = iC2 exp zk tk /0 zk+1 tk ib+ /+ zk+ tk ib+ /+ tk / Re+ Im+ exp t dt = I2 + iI2, zk + C2 sk gk /b0 + S k (0 ) exp(k /0 ), + + I2 C2 R12 R22 R12 R22, I2 C2 R12 R22 + R12 R22, Re+ Im+ Re Re Im Im Re Im Im Re (4.47) R12 = b+ /+ / tk 1/+ 1/0 + b+ /+ 2 Re, R12 = tk 1/+ 1/0 / tk 1/+ 1/0 + b+ /+ 2 Im, R22 = exp (k /0 ) cos b+ zk /+ exp k /+, Re R22 = sin b+ zk /+ exp k /+.

Im (4.48) Если 0 = +, b0 = 0, b+ = 0, то R12 = 1/ tk 1/+ 1/ Re Im R12 = 0,, R22 и R22 вычисляются, как в случае b+ = b0, по формулам (4.48).

Re Im Если 0 = +, но b0 = b+ или b0 = b+ при b+ = 0, то R12 = + /b+, Re Im R12 = 0, R22 и R22 вычисляются, как в случае + = 0, по формулам (4.48).

Re Im Если 0 = +, b+ = 0, b0 = 0, то R22 = exp k /+.

R12 = zk, Re Im Im Re R12 = R22 = 0, Третье слагаемое при 0 = + zk+ tk ib+ (zk+1 t) (tk ib0 ) (t zk ) + + exp exp I3 = iC3 dt = + zk Re+ Im+ = I3 + iI3, + C3 sk gk /b0 + S k (0 ) exp (k /0 ), + + I3 C3 R13 R23 R13 R23, C3 R13 R23 + R13 R23, Re+ Im+ Re Re Im Im Re Im Im Re (4.49) I R13 = R11, Re Im Im Re Re Re Im Im (4.50) R13 = R11, R23 = R21, R23 = R21.

Если 0 = +, b0 = b+ и b0 = 0 или 0 = + и b0 = b+ = 0, то расчеты + проводятся, как в случае интеграла I1.

284 Глава 4. Комплексное уравнение переноса Четвертое слагаемое при 0 = + zk+ tk ib+ (zk+1 t) tk (t zk ) + + exp exp I4 = C4 dt = + zk Re+ Im+ = I4 + iI4, + C4 gk /+ S k (0 ) exp (k /0 ), + + I4 C4 R14 R24 R14 R24, C4 R14 R24 + R14 R24, Re+ Im+ Re Re Im Im Re Im Im Re (4.51) I R14 = R12, Re Im Im Re Re Re Im Im (4.52) R14 = R12, R24 = R22, R24 = R22.

+ Особенности выделяются, как в случае интеграла I2.

Пятое слагаемое не зависит от 0 и b0 :

zk+ tk ib+ (zk+1 t) gk + Q+ (t) exp Re+ Im+ I5 =+ dt = I5 + iI5.

+ zk При кусочно-постоянной аппроксимации функции Q± (t) по t на интервале [zk, zk+1 ]:

Q± (t) = Q± + Q± / k+1 k получаем следующие выражения:

gk Q+ + Q+ R15 R25 R15 R25 /2, Re+ Re Re Im Im (4.53) I5 k+1 k gk Q+ + Q+ Im+ Re Im Im Re I5 R15 R25 + R15 R25 /2, k+1 k R15 = tk / tk + b+ R15 = b+ / tk + b+ 2 2 Re Im (4.54),, R25 = 1 cos b+ zk /+ exp k /+, Re R25 = sin b+ zk /+ exp k /+ ).

Im Случай b+ = 0 особенностей не дает, коэффициенты упрощаются:

R25 = 1 exp k /+.

Re Im Im Re R15 = 1/tk, R15 = R25 = 0, При линейной аппроксимации по t на [zk, zk+1 ]:

Q± (t) = [(t zk )/zk ] Q± Q± k+1 k 4.4. Расчет нерассеянной компоненты имеем = gk Q+ R15 R25 R15 R25 + Re+ Re Re Im Im I5 k gk Q+ Qk k+ R35 R45 R35 R45 + R55, Re Re Im Im Re + + zk = gk Q+ R15 R25 + R15 R25 + Im+ Re Im Im Re I5 k gk Q+ Qk k+1 Re Im Im Re Im (4.55) + R35 R45 + R35 R45 + R55, + zk tk /+ b+ /+ tk /+ + b+ /+ 2 2 Re R35 = /, R35 = 2 b+ tk / + tk /+ + b+ /+ 2 Im (4.56) /, R45 = cos b+ zk /+ exp k /+ 1, Re R45 = sin b+ zk /+ exp k /+, Im R55 = zk tk /+ / tk /+ + b+ /+ Re, R55 = zk b+ /+ / tk /+ + b+ /+ Im.

При b0 = 0 V1 (t) = V1 (zk ) gk (t zk )/0, для k = 2 K VkRe = Vk1 gk1 zk1 /0.

VkIm = 0, Re + + + В результате изменения касаются трех интегралов I1, I2 и I3. Для + вычисления I1 можно использовать формулы (4.44)–(4.46), положив b0 = 0.

Два других интеграла вычисляются заново:

а) если b0 = b+ = 0, 0 = +, то R11 = R21 = 0, Im Im R11 = 1/ tk 1/+ 1/0 R21 = exp (k /0 ) exp k /+ ;

Re Re, б) если b0 = b+ = 0, 0 = + = 0, то Im Im Re Re R11 = R21 = 0, R11 = zk, R21 = exp (k /0 ).

286 Глава 4. Комплексное уравнение переноса Если 0 = +, то zk+ tk ib+ (zk+1 t) tk (t zk ) + + I2 = C2 exp exp t dt = + zk Re+ Im+ = I2 + iI2, + C2 B12 + B22 B32 B22 B32, Re+ Re Re Re Im Im I + C2 B12 + B22 B32 + B22 B32, Im+ Im Re Im Im Re I B12 = exp (k /0 ) tk zk+1 1/+ 1/0 /d Re tk 1/+ 1/0 + b+ /+ 2 2 d2, B12 = exp (k /0 ) zk+1 b+ /+ /d 2tk b+ /+ 1/+ 1/0 /d2, Im d tk 1/+ 1/0 + b+ /+ 2 2, B22 = cos b+ zk /+ exp k /+, Re B22 = sin b+ zk /+ exp k /+, Im B32 = zk gk 1/+ 1/0 /d tk 1/+ 1/0 + b+ /+ 2 2 d2, Re B32 = zk b+ /+ /d 2tk 1/+ 1/0 b+ /+ /d2.

Im Если 0 = +, b0 = 0, b0 = b+ = 0, то zk+ + + exp k /0 + ib+ (zk+1 t) /0 t dt, I2 C zk B12 = exp (k /0 ) / b+ / Re, B12 = zk+1 exp (k /0 ) / b+ /0, Re Re Im коэффициенты B22 и B22 можно рассчитывать по формулам (4.48), когда 0 = +, B32 = 1/ b+ /0, B32 = zk / b+ /0.

Re Im Если 0 = +, b0 = b+ = 0, то + = C2 exp (k /0 ) zk+1 zk /2, 2 Re+ Im+ = I2 I 4.4. Расчет нерассеянной компоненты или можно воспользоваться формулой (4.47), положив B12 = exp (k /0 ) (zk zk+1 ) zk /2, Re Im Re Im Re Im B12 = B22 = B22 = B32 = B32 = 0.

Если + = 0, b+ = b0, то zk+ + + exp tk (t zk ) /0 tk ib+ (zk+1 t)/+ dt, I3 C3 zk + zk tk gk S k (0 ) exp (k /0 ) / 0 +, C3 + + = C3 B13 B23 B13 B23, Re+ Re Re Im Im Im+ Re Im Im Re I3 I3 = C3 B13 B23 + B13 B23, B13 = R14 = R12, B13 = R14 = R12, Re Re Im Im Im Re Re Re Re Im Im Im B23 = R24 = R22, B23 = R24 = R22.

Если + = 0, b0 = b+ = 0, то + + I3 = C3 exp (k /0 ) zk.

Остановимся на величинах b0 и b+. Обычно b0 = 0, если 0 = 0, sin 0 = 0, т. е. при нормальном падении внешнего потока, или если px = 0, т. е. уравнение переноса обычное. Величина b+ = 0, если sin + = 0 (такое значение + не входит в узлы гауссовой сетки), или sin + = 0, но px = 0 и py = 0 (в этом случае b+ = b0 = 0 и решается действительное уравнение переноса), или px = 0, sin = 0, py = 0 (узел = 0 всегда входит в азимутальную разностную сеть), или px = 0, py = 0, cos = 0 (узел = 90 всегда входит в разностную сеть по азимуту), или px = 0, py = 0, но px cos +py sin = 0, т. е.

tg = px /py (такие узлы могут входить в разностную сеть по азимуту).

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

При интегрировании уравнения (4.36) по характеристике s вы числения упрощаются, так как в правой части будут отсутствовать члены, содержащие ПЧХ нерассеянного излучения: V1 (z) 0 для s. Источник в правой части Q (t) M E P (t) g(t) Q (t) + c(t), и граничное условие 0 (zK = H) = q dV1 (H), где при z = H значения V1 (H) = VK = V1Re (H) + iV1Im (H) вычисляются по формуле (4.45), если b0 = 0. Если b0 = 0, то из формулы (4.42) находим z = g (H)/0, g (z) V1Re (H) V1Im (H) = 0, V1 (H) = g(t) dt.

288 Глава 4. Комплексное уравнение переноса Если b0 = 0, то значение V1 (H) вычисляется аналитически:

K (gk /b0 ) {sin [b0 (H zk+1 ) /0 ] sin [b0 (H zk ) /0 ]}, V1Re (H) = k= K (gk /b0 ) {cos [b0 (H zk ) /0 ] cos [b0 (H zk+1 ) /0 ]}.

V1Im (H) = k= Граничное значение 0 = q dV1 (H), т.е 0Re 0Im = q dV1Re (H), = q dV1Im (H).

K K K Интеграл zk+ gk Q (t) exp tk ib (t zk )/| | dt = I Re Im I1 + iI | | zk вычисляется с учетом способа аппроксимации функции Q (t) на [zk, zk+1 ].

При кусочно-постоянной аппроксимации = gk Q + Q R1 R2 R1 R2 /2, Re Re Re Im Im I1 k+1 k = gk Q + Q Im Re Im Im Re (4.57) I1 R1 R2 + R1 R2 /2, k+1 k R1 = tk / tk + b R1 = b / tk + b 2 2 Re Im,, R2 = 1 cos b zk /| | exp k /| |, Re R2 = sin b zk /| | exp k /| |.

Im (4.58) При линейной аппроксимации Q (t) = gk Q R1 R2 R1 R2 + Re Re Re Im Im I1 k + gk Q Q R26 R36 R26 R36 + R16 /zk, Re Re Im Im Re k+1 k = gk Q R1 R2 + R1 R2 + Im+ Re Im Im Re I1 k + gk Q Q Re Im Im Re Im (4.59) R26 R36 + R26 R36 + R16 /zk, k+1 k R16 = | | tk b 2 /a2, Re R16 = 2| |tk b /a2, a tk + b Im (4.60), R26 = cos b zk /| | exp k /| |, Re R26 = sin b zk /| | exp k /| |, Im R36 = k /a | | tk (b )2 /a2, Re R36 = b zk /a 2| |tk b /a2.

Im 4.5. Азимутальная симметрия ПЧХ, расчет интегралов столкновений Второй интеграл zk+ gk c(t) exp tk ib (t zk ) /| | dt = I I2 Re Im + iI | | zk вычисляем аналитически:

= gk S k (0 ) exp (k /0 ) R3 R4 R3 R4 /| |, Re Re Re Im Im I = gk S k (0 ) exp (k /0 ) R3 R4 + R3 R4 /| |, Im Re Im Im Re (4.61) I R3 = tk 1/| | + 1/0 / tk 1/| | + 1/0 + b /| | 2 Re, R3 = b / | | tk 1/| | + 1/0 + b /| | 2 Im, R4 = 1 cos b zk /| | exp k / 1/| | + 1/ Re, R4 = sin b zk /| | exp k / 1/| | + 1/ Im (4.62).

На первой итерации для расчета нерассеянной составляющей ПЧХ W10 (z, p, s) имеем краевую задачу Lz W 0 = 0, W0 W = 0, (4.63) = fH, 0 H 0Re+ 0Im+ которая решается точно:Wkij = 0, Wkij = 0, k = 1K;

для K = (K1) Wkij = fH exp [ (zk ) (H)] /| | cos b (H zk )/| |, 0Re (4.64) j ji j = fH exp [ (zk ) (H)] /| | sin b (H zk )/| |.

0Im Wkij j ji j § 4.5. Азимутальная симметрия ПЧХ, расчет интегралов столкновений При расчетах однопараметрических зависимостей ПЧХ либо от px, либо от py учет азимутальной симметрии или антисимметрии компонент позволяет существенно сократить объем вычислений. Анализ симметрии решения по азимуту проводится на основе метода итераций и формулы интегрирования по характеристике с детальным рассмотрением симметрии всех компонент краевой задачи: источников, коэффициентов, интегральных операторов урав нения.

При px = 0, py = 0 имеет место азимутальная симметрия W10 (z,, ) = W10 (z,, ) = W10 (z,, 2 ).

290 Глава 4. Комплексное уравнение переноса В этом случае для сокращения объема вычислений разностная сеть по азимуту вводится только для [0, ] и для расчета интеграла столкновений B(z,, ) = s (z) dex[ 2ex] W10 (z,, ) [(z,,, cos( )) + (4.65) + (z,,, cos( + ))] d используется квадратурная формула вида I J a+ Re+ [(zk, cos 1 ) + (zk, cos 2 )] + Re+ Bkj i = sk bi j kji i=1 j= J a Re [(zk, cos 1 ) + (zk, cos 2 )] +, j kji j= ± W10 (zk, ±, i ), (4.66) j kji где cos 1 = + sin sin cos( ), cos 2 = + sin sin cos( + ), a± – веса квадратуры по на интервалах + [0, 1] и + [1, 0) j– соответственно, bi – веса квадратуры по азимуту на интервале [0, ].

– Узлы разностной сети по и совпадают с узлами квадратур. Обычно по используются квадратуры Гаусса, Лобатто или специальные квадратуры для быстроосциллирующих функций, а по азимуту – формула трапеции или – квадратура Гаусса. Аналогичные выражения имеют место для компонент Re Im± Bkji, Bkji.

Если px = 0, py = 0, то действительная компонента симметрична:

W10 (z,, ) = W10 (z,, ) = W10 (z,, 2 ), Re Re Re а мнимая компонента антисимметрична по азимуту:

W10 (z,, ) = W10 (z,, ) = W10 (z,, 2 ), Im Im Im и в интегралах столкновений можно перейти к [0, ]:

Re Re ( ) [(cos 1 ) + (cos 2 )] d, B () = s (z) d Im ( ) [(cos 1 ) (cos 2 )] d.

B Im() = s (z) d Для расчета используется квадратурная формула типа (4.66).

Re Im Если px = 0, py = 0, то в каждой из компонент W10, W10 можно выделить четную и нечетную составляющие по азимуту и использовать разностную сеть 4.5. Азимутальная симметрия ПЧХ, расчет интегралов столкновений по азимуту на интервале [0, ]. Соответствующее разделение произойдет с интегралами столкновений, которые сведутся к двум выражениям типа (4.66). Но такая процедура с переходом на интервал [0, ] не имеет преимущества перед алгоритмом с заданием [0, 2], так как в обоих случаях требования к памяти ЭВМ совпадают.

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

Можно показать, что имеет место симметрия W10 (b) = W10 (b), Re Re Im Im W10 (b) = W10 (b), где b (p, s ). Это свойство позволяет сократить объем вычислений интегралов на единичной сфере или полусфере. В частности, установлено, что функция c(p) RW10 (p) при ламбертовом отражении принимает только действительные значения.

У функции W01 в общем случае px = 0, py = 0 симметрия по азимуту отсутствует уже в нерассеянной компоненте 0.

В случае px = 0, py = 0 имеем W01 (z, px,, ) = W01 (z, px,, ) = W01 (z, px,, 2 ).

При px = 0, py = W01 (z, py,, ) = W01 (z, py,, ) = W01 (z, py,, 2 ), Re Re Re W01 (z, py,, ) = W01 (z, py,, ) = W01 (z, py,, 2 ).

Im Im Im При сильной анизотропии рассеяния индикатрису представляем в виде суперпозиции (z, cos ) = c(z)(1 cos ) + (z) (z, cos ).

Интеграл столкновений на первой итерации рассчитываем по формуле S0 = s (z)c(z)0 + s (z)(z)S0, в которой квадратура берется со сглаженной индикатрисой рассеяния. Но уже со второй итерации ИМХ применяем к преобразованному уравнению (N ) + (p, z)(N ) = S(N 1) + S0, z где S (z)s (z) ds, (z, p) t (z) i(p, s ), t (z) t (z) s (z)c(z).

292 Глава 4. Комплексное уравнение переноса 4.5.1. Уменьшение размерности задачи для ПЧХ слоя с ламбертовой границей. Из соображений симметрии (среда и -источник образуют ци линдрически симметричную систему) для функции влияния горизонтально однородной атмосферы (z, x, y,, ) выполняется равенство (z, x, y,, ) = (z, x, y,, 0), (4.67) x = x cos + y sin, y = x sin + y cos.

Вводя ПЧХ W (z, px, py,, ) = (z, x, y,, ) exp i(x px + y py ) dx dy, сделаем замену переменных (4.67). Якобиан преобразований равен 1 и x = x cos y sin, y = x sin + y cos. В результате перехода к новым переменным py = px sin + py cos px = px cos + py sin, получим W (z, px, py,, ) = W (z, px, py,, 0), т. е. ПЧХ также инвариантна относительно поворота на угол. Из этого обстоятельства вытекает возможность уменьшить размерность задачи для ПЧХ: можно рассчитывать ПЧХ только для px = 0, py = 0 и всех значений азимута [0, 2], а затем пересчитывать на любые пары значений px, py.

Задача заключается в определении по значениям px, py, таких значений px, 1, чтобы W (px, py, ) = W (px, 0, 1 ).

Переходя к нулевому азимуту, получаем соотношения px cos py sin = px cos 1, px sin + py cos = px sin 1, или px = px cos(1 ), py = px sin(1 ), откуда находим 2 1/ tg(1 ) = py /px.

px = px + py, В частности, при px = 0 tg 1 = ctg, 1 = + /2. По существу это означает переход от декартовых координат, развернутых на угол, к полярным координатам. Таким образом, расчет вариаций яркости в двумерном случае сводится к интегрированию по px от 0 до и интегрированию по 1 от 4.5. Азимутальная симметрия ПЧХ, расчет интегралов столкновений до 2 (якобиан преобразования равен px ):

(px, py )q (px, py ) exp i(xpx + ypy ) dpx dpy = 10 (x, y) = (2) 1 = px(px, 1 + )q (px, 1 + ) exp [i (x cos 1 + y sin 1 )] dpx d1.

(2) 4.5.2. ПЧХ изотропно рассеивающего слоя. При фиксированной про странственной частоте p интеграл столкновений в задаче с изотропной индикатрисой рассеяния + [t (z) i(p, s )] = s (z)a(z, p), = 0, =0 (4.68) z 0 H является функцией высоты и зависит от параметра p:

2 a(z, p) = (z, p,, ) d d.

0 Интегрирование по характеристике приводит к следующим выражениям:

z + s (t)a(t, p) exp [ (t) (z) + i(z t)(p, s )] /+ dt, =+ H s (t)a(t, p) exp [ (z) (t) + i(t z)(p, s )] /| | dt + = | | z + exp [ (z) (H) + i(H z)(p, s )] /| |. (4.69) + по азимуту [0, 2] и [0, 1], Проинтегрировав выражение для а выражение для – по [0, 2] и [1, 0] и сложив эти интегралы, – получим интегральное уравнение для определения a(z) (параметр p пока опускаем):

z H (4.70) a(z) = K1 (z, t)a(t) dt + K2 (z, t)a(t) dt + f (z), z где при p = 2 (t) K1 (z, t) = s exp {[ (t) (z) + i(z t)(p, s )] /} / d d = s (t) exp {[ (t) (z)] /} J0 (z t) 1 2 p2 + p2 / / d. (4.71) = x y 294 Глава 4. Комплексное уравнение переноса Аналогично (t) K2 (z, t) = s exp {[ (z) (t)] /||} J0 (t z) 1 2 p2 + p2 /|| /|| d, (4.72) x y exp {[ (z) (H)] /||} J0 (H z) 1 2 p2 + p2 /|| d.

f (z) = x y В случае однородного слоя (z) = t z и уравнение (4.70) превращается в уравнение типа свертки. Если при этом p = 0, то приходим, естественно, к обычному интегральному уравнению Пайерлса для однородного, изотропно рассеивающего плоского слоя с единичным изотропным источником на границе:

H a(z) = s a(t)E1 (t |z t|) dt + E [t (H z)].

Нетрудно определить предельные значения + и на границах слоя при p = 0:

lim + (H) = (s /t ) a(H) при + 0, lim (0) = (s /t ) a(0) при 0.

Оценим асимптотику ПЧХ изотропно рассеивающего слоя. Для этого рассмотрим асимптотическое поведение f (z) при |p|. Произведя замену переменной t = 1/||, получим exp {[ (z) (H)] t} J0 (H z) t2 1 p2 + p2 /t2 dt.

f (z) = x y Покажем, что существует константа C, такая, что C t2 1 dt, (4.73) f (z) exp(1 t)J где 1 (H) (z), (H z) p2 + p2.

x y Действительно, условие (4.73) равносильно неравенству exp (1 t) C 1/t2 J0 t2 1 dt 0. (4.74) Принимая во внимание свойства функции Бесселя J0, для выполнения (4.74) достаточно показать, что существует такое C, что функция v(t) = 4.6. Алгоритм расчета пространственно-частотной характеристики exp(1 t) C 1/t2 монотонно убывает при t 1. Производная v (t) = exp(1 t) 1 C 1/t2 + 2/t отрицательна при C (1 + 2) /1. В этом случае справедлива оценка C exp 12 + 2 / 12 + f (z, p) и асимптотика ПЧХ имеет порядок O exp 12 + 2 / 12 + 2, поскольку SW1 = s (z)f (z, p).

Интегральные уравнения второго рода для a(z) (4.70) решаются численно с помощью стандартных алгоритмов, в частности, методом квадратур. Сами ПЧХ находим с помощью численного интегрирования исходного уравнения (4.68) по характеристике (4.69).

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

Ядрами этих функционалов – оптических передаточных операторов – – – являются функции влияния (ФВ) краевых задач теории переноса или пространственно-частотные характеристики (ПЧХ) – фурье-образ ФВ по гори – зонтальным координатам. Другими словами, решения общей и первой краевых задач теории переноса находятся с помощью фундаментальных решений, или ФВ, которые определяются методом преобразования Фурье. Этот подход, называемый методом ФВ и ПЧХ, позволяет решать прямые задачи (на его основе разрабатываются радиационные блоки в разных прикладных моделях), а также создавать новые методики решения обратных задач в проблемах дистанционного зондирования.

Метод ФВ и ПЧХ развит на случай ортотропной подстилающей поверхно сти, анизотропной отражающей и пропускающей границы раздела двух сред (например, атмосфера–океан), а также с учетом поляризации излучения.

296 Глава 4. Комплексное уравнение переноса Метод ФВ и ПЧХ позволяет редуцировать трехмерную по пространству, пятимерную по фазовому объему краевую задачу для кинетического инте гродифференциального уравнения к параметрическому набору одномерных по пространству, трехмерных по фазовому объему краевых задач для непре рывной, бесконечно дифференцируемой и убывающей по действительному параметру комплекснозначной функции, т. е. параметрических одномерных комплексных уравнений переноса, различающихся источниками в зависимости от вида «возмущения».

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

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

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

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

радиус-вектор r = (x, y, z). На верхнюю границу слоя z = 0, неограниченного в горизонтальном направлении ( x, y ) и конечного по высоте H), под углом 0 = arccos 0 падает внешний параллельный (0 z солнечный поток мощности S в азимутальной плоскости 0, s0 = (0, 0 ).

Предполагается, что рассеяние в акте происходит без изменения частоты.

Подстилающая поверхность на границе слоя z = H отражает падающее изнутри излучение по заданному закону. Требуется определить интенсив 4.6. Алгоритм расчета пространственно-частотной характеристики ность многократно рассеянного излучения внутри слоя, а также излучения, отраженного и пропущенного слоем.

Множество всех направлений s = (, ), = cos [1, 1], [0, 2], образует единичную сферу = +, где + и – полусферы для – направлений распространения нисходящего, пропущенного излучения ( [0, 1]) и восходящего, отраженного излучения ( [1, 0]) соответственно.

Для удобства записи граничных условий вводим множества t = {(z, r, s) : z = 0, s + }, b = {(z, r, s) : z = H, s };

r = (x, y) – проекция радиус-вектора на горизонтальную плоскость. Для – удобства и большей наглядности величины, определенные для направлений s = s+ +, помечаются меткой «+», а для s = s – меткой «».

– В предположении стационарного состояния среды и постоянства источника инсоляции поле квазимонохроматического оптического и миллиметрового излучения описывается интенсивностью излучения (r, s). Для макроскопи чески изотропной и плоскостратифицированной немультиплицирующей (без размножения) системы атмосфера–подстилающая поверхность интенсивность находится как решение стационарной общей краевой задачи для интегро дифференциального уравнения переноса с линейными операторами: оператор переноса D Dz + sin cos Dz (4.75) + sin sin, + tot (z);

x y z интеграл столкновений S sc (z) (z, s, s )(z, r, s ) ds, (4.76) где ds = d d ;

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

В качестве коэффициентов операторы содержат физические характери стики среды: вертикальные распределения суммарного коэффициента ослаб ления (экстинкции) tot (z) = sc (z) + a (z), суммарных коэффициентов рассеяния sc (z) = a,s (z) + R (z) и поглощения a (z) = a,a (z) + R,a (z) с аэрозольными (a,s (z) и a,a (z)) и молекулярными (R (z) и R,a (z)) компонентами, суммарной индикатрисы рассеяния a,s (z) (z) a (z, ) + R (4.77) (z, ) = R (), sc (z) sc (z) которая определяется через индикатрисы аэрозольного a (z, ) и рэлеевского рассеяния (1 + cos2 ), R () = 298 Глава 4. Комплексное уравнение переноса с нормировкой ds = 2 (z, cos ) d cos = 1;

(4.78) – угол рассеяния из направления s = (, ) в направление s = (, ).

– Если источник анизотропный горизонтально-неоднородный и при этом не расщепляются пространственные и угловые переменные, то функция влияния (sH ;

z, r, s) – решение краевой задачи – = f (sH ;

r, s) K = 0, = 0, (4.79) t b с параметром sH и источником f (sH ;

r, s) = (r )(s sH ).

Если источник изотропный горизонтально-неоднородный, то функция влияния (s ;

z, r, s) ds (4.80) r (z, r, s) = удовлетворяет краевой задаче K = 0, = 0, (4.81) = (r ).

t b Функция влияния (sH ;

z, r, s) описывает фактически поле излучения в слое с неотражающими границами, создаваемое за счет процессов многократ ного рассеяния стационарного лазерного луча с направлением sH, источник которого расположен на границе z = H в центре системы горизонтальных координат x, y. Функция влияния r (z, r, s) совпадает с функцией размытия точки;

ее фурье-образ по r в направлении надира, когда s = ( = 1, = 0), совпадает с частотно-контрастной характеристикой.

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


F[f (r )](p) = (4.82) f (r ) exp[i(p, r )] dr.

Первая краевая задача для трехмерного уравнения переноса (4.79) с помощью фурье-преобразования (4.82) приводится к первой краевой задаче для параметрического одномерного уравнения переноса с оператором L(p) Dz i(p, s ) S, (p, s ) = px sin cos + py sin sin.

Пространственная частота p = (px, py ) принимает только действительные значения ( px, py );

s = (sin cos, sin sin ) – проекция вектора – 4.6. Алгоритм расчета пространственно-частотной характеристики направления s на горизонтальную плоскость. Пространственно-частотная характеристика (sH ;

z, p, s) – решение комплексного уравнения переноса – = g (sH ;

p, s) L(p) = 0, = 0, (4.83) t b с источником g (sH ;

p, s) = F[f (sH ;

r, s)] = (s sH ), так как F[(r ] = 1. ФВ и ПЧХ связаны прямым и обратным фурье преобразованиями:

= F 1 [].

= F[], ФВ r отвечает ПЧХ r (z, p, s) = F[r ], которая удовлетворяет краевой задаче L(p)r = 0, r = 0, r = 1 (4.84) t b с единичным изотропным источником.

Численное решение краевой задачи (4.79) или (4.81) для ФВ воз можно только алгоритмами метода Монте-Карло. Для применения конечно разностных методов необходимы специальные меры для ограничения размера слоя в горизонтальной плоскости. Наиболее эффективным и естественным оказывается подход, основанный на определении ФВ методом фурье-преобра зования. При этом от задачи для трехмерного, неограниченного по горизон тальным координатам слоя осуществляется переход к задаче для одномерного, конечного по высоте слоя. Одновременно происходит факторизация решения:

вместо координаты r появляется параметр p и задача для ПЧХ решается при фиксированных значениях этого параметра.

Краевая задача для ПЧХ (4.83) или (4.84) формально совпадает с краевой задачей теории переноса излучения в одномерном плоском слое с тем лишь принципиальным различием, что искомая функция (z, p, s) и коэффициент экстинкции (p, z) tot (z)i(p, s ) являются комплекснозначными и зависят от действительного параметра p.

4.6.2. Выделение фазового множителя. В силу линейности краевой за дачи, в ПЧХ можно выделить две составляющие: = 0 + d. Прямая, нерассеянная компонента 0 от источника находится как решение задачи M (p)0 = 0, 0 = g (s0 ;

s), 0 = g (sH ;

s) (4.85) t b интегрированием по характеристикам дифференциального оператора M (p) Dz i(p, s ) в явном виде:

iz(p, s+ ) (z) 0 (z, p, s+ ) = g (s0 ;

s+ ) exp (4.86) exp, + + i(H z)(p, s ) (H) (z) 0 (z, p, s ) = g (sH ;

s ) exp (4.87) exp.

| | | | 300 Глава 4. Комплексное уравнение переноса Диффузная, многократно рассеянная компонента d (z, p, s) удовлетворяет краевой задаче L(p)d = B 0, d = 0, d = 0 (4.88) t b S0, содержащим B 0 (z, p, s) с внутренним, распределенным источником осциллирующие множители, из-за которых возникают осцилляции в решении задачи (4.88). Такие осцилляции осложняют расчет интеграла столкновений методом квадратур.

Принимая во внимание выражения (4.86) и (4.87), выделим фазовые множители в решении задачи (4.83):

iz(p, s+ ) (z, p, s) = + (z, p, s+ ) exp s = s+ + ;

(4.89), i(H z)(p, s ) (z, p, s) = (z, p, s ) exp s = s.

(4.90), || Вместо задачи (4.83) для осциллирующей функции будем численно решать краевую задачу для сглаженных функций + и :

Dz + = B +, s = s+ + ;

(4.91) Dz = B, s = s ;

(4.92) + = 0, = g (sH ;

s) (4.93) t b с функциями источников iz(p, s+ ) iz(p, s+ ) B + exp S + + exp + + + i(H z)(p, s ) + S exp = e1 S + + e3 + S e4 ;

(4.94) | | i(H z)(p, s ) iz(p, s+ ) B exp S + + exp + | | + i(H z)(p, s ) + S exp = e2 S + + e3 + S e4.

(4.95) | | Для краткости записи введены обозначения:

iz(p, s+ ) i(H z)(p, s ) e1 exp e2 exp ;

;

+ | | iz(p, s+ ) i(H z)(p, s ) e3 exp e4 exp ;

.

+ | | Интегралы столкновений определены на полусферах:

S + sc S sc ds, ds.

+ 4.6. Алгоритм расчета пространственно-частотной характеристики Дифференциальный оператор переноса в уравнениях (4.91)–(4.92) принял форму, совпадающую с обычным уравнением переноса.

4.6.3. Метод подобия. Среди приближенных методов решения одномерного уравнения переноса достаточно распространен метод подобия, когда в качестве пространственной переменной берется оптическая толщина. В этом случае от уравнения или (4.96) + (, s) = ( ) ds, D = S, где альбедо акта рассеяния = sc /tot, при выделении сильной анизотропии рассеяния в виде (cos ) = c(1 cos ) + (1 c) (cos ) (4.97) для вертикально однородного слоя можно перейти к уравнению (4.98) + (, s) = ( ) ds, где новая оптическая толщина (1 c) = (1 c) ;

(4.99) = (1 c) – новое значение альбедо акта рассеяния. С помощью соотношений подо – бия (4.99) вместо уравнения (4.96) численно решается уравнение (4.98) со сглаженной индикатрисой рассеяния (cos ). В наиболее широко ис пользуемом приближенном методе -Эддингтона в качестве (cos ) берется двухчленное разложение по полиномам Лежандра.

В уравнениях (4.91)–(4.92) для определения ПЧХ дифференциальный опе ратор можно записать через оптическую толщину в качестве пространственной переменной:

iz(p, s+ ) iz(p, s+ ) D + = exp S + exp + + + + i(H z)(p, s ) + S exp (4.100) ;

| | i(H z)(p, s ) iz(p, s+ ) D = exp S + exp + + | | + i(H z)(p, s ) + S exp (4.101), | | с интегральными операторами + S ( ) S ( ) ds, ds.

+ 302 Глава 4. Комплексное уравнение переноса Однако, важно отметить, что наряду с оптической толщиной задачи (4.100)– (4.101) содержат в фазовых множителях пространственную переменную z.

Если индикатриса рассеяния представляется в виде (4.97), то вме сто (4.100)–(4.101) можно получить преобразованные уравнения, которые + по форме совпадут с (4.100)–(4.101), если операторы D, S, S заменить соответственно на операторы + D c + (1 c), S c (1 c) ds, + + S c (1 c) ds, S c = S c + S c.

Как и в случае обычного уравнения переноса, для однородного слоя можно перейти к и согласно соотношениям (4.99). Однако, метод подобия для решения комплексного уравнения переноса не применим, так как в фазовых множителях сохраняется зависимость от координаты z.

Такой же вывод вытекает из рассмотрения комплексного уравнения переноса (4.83), в котором не выделен фазовый множитель. Действительно, если в уравнении (4.83) перейти к оптической толщине:

i(p, s ) D = S tot ( ) и воспользоваться представлением индикатрисы рассеяния (4.97), то получим i(p, s ) D c (4.102) = S c.

tot ( ) При введении и согласно (4.99) в левой части преобразованного уравнения (4.102) сохраняется слагаемое с комплексным коэффициентом, равным i(p, s )/[tot (1 c)], что не позволяет использовать метод подобия.

Уравнения (4.100)–(4.101) совпадают по структуре с уравнениями (4.91)– (4.92). Оператор переноса D можно получить из оператора Dz, если положить z = и tot (z) = 1, а переход от S к S осуществляется путем замены sc (z) на ( ). Далее будем излагать алгоритм численного решения краевой задачи (4.91)–(4.92).

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

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

a (z, ) = c(z)(1 cos ) + (z)a (z, ). (4.103) Полная исходная аэрозольная индикатриса рассеяния a, bx (z, ) перенорми руется по формуле (z, ) a (z, ) = a, bx (4.104) ;

Aa, bx (z) = a, bx (z, ) d cos.

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

a, bx (z, ) (4.105) a (z, ) = ;

Aa, bx (z) = a, bx (z, ) d cos.

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

f (z) A (z) (1 cos ) + a, bx (4.106) a (z, ) = a (z, ).

2Aa, bx (z) Aa, bx (z) Откуда вытекает, что f (z) = Aa, bx (z) Aa, bx (z) 0, причем f (z) = 0, если «сглаженная» аэрозольная индикатриса рассеяния совпадает полностью с исходной аэрозольной индикатрисой рассеяния. Со поставляя выражения (4.103) и (4.106), получаем Aa, bx (z) f (z) [1 (z)].

(z) =, c(z) = = 2Aa, bx (z) Aa, bx (z) Суммарная аэрозольно-молекулярная индикатриса рассеяния (4.77) с выделением -анизотропии представляется в виде, аналогичном (4.103):

(z, ) = v(z)(1 cos ) + (z, ) (4.107) с коэффициентом a,s (z) v(z) = c(z) a,s (z) + R (z) 304 Глава 4. Комплексное уравнение переноса и сглаженной составляющей a,s (z) R (z) (z, ) = (z)a (z, ) + R ().

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


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

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

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

4.6.5. Последовательность вычислений с учетом -анизотропии рассея ния и источника. Опишем кратко последовательность вычислений ПЧХ как решения краевой задачи (4.91)–(4.93) с двумя типами источников, отвечающим изотропным, когда g (sH ;

s) = 1, и анизотропным, когда g (sH ;

s) = (sH s), граничным условиям.

Выделим прямое излучение + = + + +, = + (4.109) пр пр d d как решение задач Коши Dz + = 0, + = 0, s = s+ + ;

(4.110) пр пр t Dz = 0, s = s, = g (sH ;

s), (4.111) пр пр b которые интегрированием по характеристикам разрешаются в явном виде:

+ (z, p, s) = 0, s = s+ + ;

(4.112) пр (H) (z) (z, p, s) = g (sH ;

s) exp s = s. (4.113), пр || 4.6. Алгоритм расчета пространственно-частотной характеристики С учетом представления индикатрисы рассеяния (4.108) и соотношения (1 cos ) = (s s ) правые части (4.94)–(4.95) можно условно записать в виде следующих выражений:

+ B + B + + Bс ;

+ (4.114) B B (4.115) + Bс.

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

+ B + a, s (z)c(z)+ (z, p, s), s = s+ + ;

(4.116) B a, s (z)c(z) (z, p, s), s = s, (4.117) + и со «сглаженной» индикатрисой, причем выражение для Bс совпадает – с (21), если в интегралах столкновений S + и S с (4.94), а для Bс – + заменить (z, ) на (z, ) и ввести обозначения Sс и Sс соответственно.

В задаче для многократно рассеянной, диффузной компоненты ПЧХ Dz + = B + d + B + пр, (4.118) d Dz = B d + B пр, (4.119) d + = d = 0, (4.120) db t источники определяются через решения (4.112)–(4.113) с учетом представле ний (4.94)–(4.95):

F0 (z, p, s) B + пр = e1 S e4 ;

+ (4.121) пр F0 (z, p, s) B пр = e2 S e4.

(4.122) пр При анизотропном граничном условии в задаче (4.91)–(4.93) для источ ников (4.121)–(4.122) получаются явные выражения:

+ F0 = e1 eH e sc (z)(z, cos H );

(4.123) F0 = e2 eH e sc (z)(z, cos H ), (4.124) где для краткости записи введены обозначения:

i(H z)(p, sH ) (H) (z) eH exp e exp ;

;

|H | |H | cos H = H + sin sin H cos( H ).

Решение задачи (4.118)–(4.120) методом последовательных приближений – – методом простых итераций по кратности рассеяния есть сумма ряда Неймана:

N N + = + ;

=, (4.125) n n d d n=1 n= 306 Глава 4. Комплексное уравнение переноса компоненты которого определяются из системы рекуррентных задач:

+ Dz + = Fn1, Dz = Fn1, + = 0, =0 (4.126) n n n n t b с источниками + Fn1 B + n1, Fn1 B n1. (4.127) Учитывая специальным образом -анизотропию рассеяния и источника, предлагаем последовательность вычислений, отличную от простых итера ций (4.126).

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

N d = 1 + d1, где d1 = n.

n= Этап 1. Полный вклад однократного рассеяния находится из задачи Dz + = F0, + Dz = F0, + = 0, =0 (4.128) 1 1 1 t b с источниками (4.123)–(4.124), содержащими полную аэрозольно-рэлеевскую индикатрису рассеяния. С учетом -анизотропии рассеяния (4.107) функции источников (4.123)–(4.124) преобразуются к виду:

+ + F0 = F0, (1 cos H ) + F0, c, (4.129) F0 = F0, c ;

где гладкие составляющие + F0, c e1 eH e sc (z) (z, cos H );

(4.130) e2 eH e sc (z) (z, cos );

H (4.131) F0, c коэффициент в сингулярной компоненте F0, a, s (z)c(z)e. (4.132) Вследствие линейности краевой задачи (4.128) с учетом представления источника (4.129) можно разложить 1 на две компоненты:

+ = + c ;

= 1, (1 cos H ) + c. (4.133) 1 1, 1 1, Этап 2. Вклад однократного рассеяния со сглаженной индикатрисой рассеяния – гладкая составляющая ПЧХ – находится из задачи – – Dz + c = F0, c, + Dz c = F0, c, + c = 0, c =0 (4.134) 1, 1, 1, 1, t b с гладкими функциями источников (4.130)–(4.131).

Этап 3. Сингулярная составляющая однократного рассеяния рассчиты вается с помощью решения задачи s = sH, = 0, (4.135) DH 1, = F0,, 1, b 4.6. Алгоритм расчета пространственно-частотной характеристики где оператор DH совпадает с оператором Dz при = H. В задаче для расчета d1 :

Dz + = B + d1 + B + 1, Dz = B d1 + B 1, d1 d + d1 = 0, = d t b проведем преобразование коэффициента экстинкции, переходя к представле ниям (4.114)–(4.115) для функций источника:

Dzc + = Bс d1 + B + 1, Dzc = Bс d1 + B 1, + d1 d + d1 = 0, = 0;

d t b оператор переноса Dzc + [tot (z) a, s (z)c(z)]. (4.136) z Выделяем вклад второго рассеяния:

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

n= Этап 4. Вклад второго рассеяния находится как решение задачи Dzc + = F1+, Dzc = F1, + = 0, =0 (4.137) 2 2 2 t b с источниками, которые определяются через вклад полного однократного рассеяния – решение задачи (4.128) по формулам (4.94)–(4.95):

– F1+ B + 1 = e1 S + + e3 + S e4 ;

(4.138) 1 F1 B 1 = e2 S + + e3 + S e4, (4.139) 1 или с учетом представлений (4.114)–(4.115) F1+ = B + + Bс 1 = a, s c(z)+ c + + + 1 1, + e1 e4 1, sc (z) (z, cos H ) + e1 Sс + c e3 + Sс c e4 ;

+ (4.140) 1, 1, F1 = B + Bс 1 = a, s c(z) + 1 + e2 e4 1, sc (z) (z, cos H ) + e2 Sс + c e3 + Sс c e4.

+ (4.141) 1, 1, Учитывая -анизотропию в индикатрисе рассеяния (4.107) и в однократном приближении (4.133) в источниках (4.138)–(4.139) можно выделить сингу лярную компоненту:

F1+ = F1, c ;

+ F1 = F1, (1 cos H ) + F1, c.

(4.142) + Гладкая функция F1,c совпадает с (4.140), гладкая составляющая F1, c a, s c(z) c + e2 e4 1, sc (z) (z, cos H ) + 1, + e2 Sс + c e3 + Sс e4 ;

+ (4.143) 1, 1,c 308 Глава 4. Комплексное уравнение переноса коэффициент в сингулярной компоненте F1, a, s (z)c(z)1,. (4.144) Исходя из представления источников (4.142), второе приближение можно записать в виде суперпозиции + = + c ;

= (1 cos H ) + c. (4.145) 2 2, 2 2, 2, Этап 5. Гладкая составляющая вклада второго рассеяния рассчитывается как решение задачи Dzc + c = F1, c, + Dzc c = F1, c, + c = 0, c =0 (4.146) 2, 2, 2, 2, t b с гладкими функциями источников (4.140) и (4.143).

Этап 6. Сингулярная компонента во втором рассеянии находится с помощью решения задачи s = sH ;

= 0, (4.147) DHc 2, = F1,, 2, b оператор DHc совпадает с оператором Dzc при = H.

Разложение (4.145) используется для расчета источника – сглаженного – интеграла столкновений + + Bс 2 = Bс 2, c + e1 e4 2, sc (z) (z, cos H ), (4.148) H (4.149) Bс 2 = Bс 2, c + e2 e4 2, sc (z) (z, cos ) в задаче, описывающей вклад диффузного излучения, начиная с третьей кратности рассеяния:

Dzc + = Bс d2 + Bс 2, + + d Dzc = Bс d2 + Bс 2, d + d2 = 0, = 0. (4.150) d t b Этап 7. Расчет вклада третьего рассеяния как решение задачи Dzc + = Bс 2, Dzc = Bс 2, + = 0, + =0 (4.151) 3 3 3 t b с гладкими функциями источников (4.148)–(4.149).

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

Dzc +N = Bс N 1 + Bс 2, + + d2 d Dzc N = Bс N 1 + Bс 2, d2 d +N = 0, N = 0. (4.152) d2 d t b 4.6. Алгоритм расчета пространственно-частотной характеристики При изотропном граничном условии в задаче (4.91)–(4.93) выраже ния (4.112)–(4.113) не содержат сингулярностей. С учетом -анизотропии рассеяния (4.107) получаются гладкие функции источников:

F0 B + пр = e1 Sс e4 ;

+ (4.153) пр F0 B пр = a, s (z)c(z)eH + e2 Sс e4.

(4.154) пр В этом случае задачу для диффузной компоненты ПЧХ (4.118)–(4.120) можно решать итерационно по схеме (n 1 – номер итерации):

– Dzc + = Bс n1 + F0, + + n Dzc = Bс n1 + F0, n + = 0, = 0, (4.155) n n t b + причем для n = 1 в правой части берутся только источники F0 и F0.

Следует отметить, что в случае изотропного источника можно не выделять гладкую компоненту вклада прямого, нерассеянного излучения и решать непо средственно задачу (4.91)–(4.93) с учетом -анизотропии рассеяния (4.107) итерационно по схеме (n 1 – номер итерации):

– Dzc + = Bс n1, + Dzc = Bс n1, + = 0, = 1. (4.156) n n n n t b 4.6.6. Разностная сеть и дискретная аппроксимация задачи. Описание физических характеристик среды, составляющих оптическую модель слоя, формализуем в соответствии с едиными требованиями к дискретным алго ритмам ИМХ.

Для учета пространственной неоднородности среды по высоте в виде плоскопараллельных слоев вводятся зоны двух типов, а именно:

1) зоны с равномерным шагом пространственной разностной сети;

их границы обозначаются через zl1, l1 = 1 L1, zl1 =1 = 0, zl1 zl1 +1, zl1 =L1 = H;

шаги разностной сети в пределах l1 -зоны zl1 = (zl1 +1 zl1 )/Nl1, где Nl1 – – количество равных интервалов по z в l1 -зоне;

размерность пространственной 1 разностной сети zk, k = 1 K, z1 = 0, zk=K = H, равна K = L=1 Nl1 + 1;

l шаг разностной сети zk = zk+1 zk = zl1, если zl1 zk zk+1 zl1 +1, т. е.

отрезок [zk, zk+1 ] лежит в пределах l1 -зоны, совпадающей с отрезком [zl1, zl1 +1 ];

2) зоны с постоянными по z индикатрисами рассеяния (z, cos );

их границы обозначаем через zl, l = 1 L, zl=1 = 0, zl zl+1, zl=L = H;

для направлений с s = s+ полагаем (z, cos ) = l (cos ), если zl z zl+1 ;

для направлений с s = s полагаем (z, cos ) = l (cos ), если zl z zl+1 ;

не нарушая общности физических постановок задачи, считаем, что границы l-зон входят в число границ l1 -зон, т. е. L1 L.

Высотный ход коэффициентов экстинкции, аэрозольного и рэлеевского рассеяния и поглощения задается послойно в виде кусочно-постоянной, 310 Глава 4. Комплексное уравнение переноса кусочно-линейной, кусочно-экспоненциальной или другой кусочно-аналитиче ской аппроксимации. Для всех точек разностной сети zk формируются массивы + коэффициентов рассеяния sc, k для s = s+, индекс которого определяется из условия zk z zk+1, и sc, k для s = s, индекс которого определяется из условия zk z zk+1. Эти массивы вводятся для учета разрывов коэффици ентов sc на границах k-интервалов в вычислениях интегралов столкновений:

при интегрировании по s + в точке zk берутся значения sc и (z, cos ) слоя, прилегающего сверху к zk, а по s – слоя, прилегающего снизу к zk.

– Для дискретной аппроксимации задачи, построения конечно-разностных схем метода интегрирования по пространственным характеристикам и квадра турных формул вводим разностную сеть = {zk, j, i }: по пространственной переменной zk, k = 1 K;

по зенитному углу + (0, 1], j = 1 J +, и j [1, 0), j = 1 J ;

по азимуту i [0, 2], i = 1 I. Сеточные функции j при фиксированном значении параметра определяем по правилу:

+ (zk, +, i ), (zk,, i ).

j j kji kji 4.6.7. Алгоритм интегрирования по характеристикам. Основы вычис лительного алгоритма интегрирования по характеристикам рассмотрим на примере краевой задачи Dz + = F +, Dz = F, + = f, (4.157) =g t b с известными комплекснозначными функциями источников F + = F Re+ + iF Im+, F = F Re + iF Im, f = f Re + if Im, g = gRe + igIm. (4.158) Характеристики одномерного дифференциального оператора Dz совпадают с направлениями распространения излучения и определяются соотношениями = const, = const. Явные формулы интегрирования по характеристикам дифференциального оператора одномерного параметрического комплексного уравнения переноса (4.157) получаются по аналогии со случаем обычного одномерного кинетического уравнения:

z z 1 F + (t, p, +, ) exp tot (u) du dt + + (z, p, +, ) = + + t z + f (s0 ;

p, +, ) exp tot (u) du = (4.159) + z (z) (t) 1 (z) F + (t, p, +, ) exp dt + f (s0 ;

p, +, ) exp = ;

+ + + 4.6. Алгоритм расчета пространственно-частотной характеристики H t 1 F (t, p,, ) exp tot (u) du dt + (z, p,, ) = | | | | z z H H 1 + g(sH ;

p,, ) exp tot (u) du = F (t, p,, ) | | | | z z (t) (z) (H) (z) dt + g(sH ;

p,, ) exp exp (4.160).

| | | | Поскольку задача (4.157) решается при фиксированных значениях пара метров p, s0, sH, то для сокращения записей в дальнейшем изложении эти параметры опускаются в перечне аргументов функций источника и ПЧХ.

Если вычисления проводить, используя комплексную арифметику, то фор мально конечно-разностный аналог выражений (4.159)–(4.160) на разностной сети = {zk, j, i } совпадает со случаем обычного уравнения переноса.

Для = + 0 на границе z = 0 (k = 1) полагаем j + = f (+, i ) = fji, (4.161) j 1ji а для k = 2 K, пользуясь свойством аддитивности экспонент, приходим к представлению zk+ + ji = + exp tot (u) du + + k+1, kji j zk zk+1 zk+ 1 F + (t, +, i ) exp + tot (u) du dt = ++ j j j zk t zk+ k+1 (t) k = + exp F + (t, +, i ) exp (4.162) + dt;

+ + + j kji j j j zk zk+1 zk+1 zk tot (u) du tot (u) du = k+1 k.

k = tot (u) du = zk 0 Аналогично для = 0 на границе z = H (k = K) полагаем j = g(, i ) = gji, (4.163) j Kji 312 Глава 4. Комплексное уравнение переноса а для k = (K 1) 1 имеет место представление zk+ ji exp = tot (u) du + kji k+1, |j | zk zk+1 t 1 (t,, i ) exp + F tot (u) du dt = j |j | |j | zk zk zk+ (t) k k = ji exp F (t,, i ) exp (4.164) + dt.

| | j k+1, |j | |j | j zk Представления (4.159) и (4.161)–(4.162) для 0, а также (4.160) и (4.163)–(4.164) для 0 взаимно однозначны на сети zk. Но вторые пред ставления в виде рекуррентных соотношений позволяют ввести существенную экономию объема вычислений и организовать последовательность расчета с переходом к следующему уровню по z от предыдущего уровня, отстоящего на расстоянии шага разностной сети zk, а не от самой границы слоя. Чтобы получить конечно-разностные схемы, необходимо ввести аппроксимацию коэффициента ослабления tot и функции источника F (z,, ) на отрезках [zk, zk+1 ]. Если функция источника F (z,, ) выражена явно, то при кусочно постоянной и некоторых других аппроксимациях коэффициента tot (z) инте гралы (4.159)–(4.160) могут быть вычислены в явном виде. Наиболее часто применяется кусочно-постоянная аппроксимация на отрезке [zk, zk+1 ]:

(4.165) tot (z) = tot, k, sc (z) = sc, k, F (z, s) = [F (zk+1, s) + F (zk, s)]/2.

4.6.8. Конечно-разностная схема расчета приближения однократного рассеяния методом характеристик. Рассматривается задача Dz + = F +, Dz = F, + = 0, = 0, (4.166) t b которая совпадает с задачей (4.128), если положить + = +, =, + F + = F0, F = F0, 1 или с задачей (4.134), если положить + = + c, = c, + F + = F0, c, F = F0, c.

1, 1, Источники в задачах (4.128) и (4.134) представлены явными выражениями, отличающимися только входящими в них индикатрисами рассеяния и соответственно. Достаточно получить решение задачи (4.166), совпадающей с (4.128), чтобы найти и решение задачи (4.134) путем замены на.

С помощью дискретного представления формул интегрирования по ха рактеристикам (4.161)–(4.162) и (4.163)–(4.164), явного выражения функций 4.6. Алгоритм расчета пространственно-частотной характеристики источника (4.123)–(4.124):

iz(p, s+ ) i(H z)(p, sH ) + F0 = exp exp + |H | (H) (z) exp sc (z)(z, cos H );

(4.167) |H | i(H z)(p, s ) i(H z)(p, sH ) F0 = exp exp | | |H | (H) (z) exp sc (z)(z, cos H ) (4.168) |H | и кусочно-постоянной аппроксимации на отрезке [zk, zk+1 ] коэффициентов уравнения переноса:

(z, cos H ) = l (cos H ) (4.169) tot (z) = tot, k, sc (z) = sc, k, получаем следующие двухточечные схемы расчета компоненты ПЧХ, обу словленной однократным рассеянием, – решения задачи (4.166) с источни – ками (4.167)–(4.168):

для k = 1 K + ji = + exp +k + k+1, kji j 1 1 + sc, k l (cos H ) tot, k + i ji + H +H + | | j + 1 1 + + H tot, k +H + ji + j + | | j (H) k + +k + i H (H zk ) ji zk exp |H | j 1 1 + izk H + ji 1 ;

(4.170) exp k ++ |H | j если = H, то для k = K 1 j k = ji exp + | | kji k+1, j 1 1 + i H ji sc, k l (cos H ) tot, k |H | | | j + 1 1 | | + H ji tot, k j | | |j | H (H) k exp + i(H zk ) H ji |H | 1 1 + izk ji H 1 ;

(4.171) exp k |H | |j | 314 Глава 4. Комплексное уравнение переноса если = H, то cos H = 1 и для k = K 1 1 и j = j H j (H) k k = ji exp + sc, k l (cos H ) exp (4.172), |H | |H | kji k+1, так как в этом случае (H) (z) F0 = sc (z)(z, cos H ) exp.

|H | Для краткости введены обозначения (p, s+ ) (p, s ) (p, sH ) + ji ji H ;

;

.

+ | | |H | j j Аналогично получаются явные выражения для расчета сингулярных компонент 1, и 2, как решений задач (4.135) и (4.147) соответственно.

4.6.9. Конечно-разностная схема интегрирования комплексного урав нения переноса c дискретным источником методом характеристик.

Рассматривается задача (4.166) с правой частью уравнений, заданной сеточной функцией источника:

F + = B + = F + (zk, +, i ) = Fkji, + j F = B = F (zk,, i ) = Fkji.

(4.173) j В комплекснозначных функциях выделяем действительные и мнимые части:

+ Re+ Im+ Re Im Fkji = Fkji + iFkji, Fkji = Fkji + iFkji, + = Re+ + iIm+, = Re + iIm. (4.174) kji kji kji kji kji kji Из дискретного представления формул интегрирования по характеристи кам (4.161)–(4.162) и (4.163)–(4.164) очевидно, что эти выражения расщеп ляются для действительных и мнимых частей:

Re+ = 0, Im+ = 0;

(4.175) 1ji 1ji для k = 1 K zk+ tot, k (zk+1 t) Re+ ji = Re+ exp +k Fji (t) exp Re+ ++ dt;

+ k+1, kji j j j zk (4.176) zk+ tot, k (zk+1 t) k Im+ji = Im+ exp Fji (t) exp Im+ + dt;

+ + + k+1, kji j j j zk (4.177) Re Im = 0, = 0;

(4.178) Kji Kji 4.6. Алгоритм расчета пространственно-частотной характеристики для k = K 1 zk+ tot, k (t zk ) k + Fji (t) exp Re Re ji exp Re = dt;

| | kji k+1, |j | |j | j zk (4.179) zk+ tot, k (t zk ) k Im = Imji exp Fji (t) exp Im + dt.

| | | | kji k+1, |j | j j zk (4.180) По аналогии со случаем обычного одномерного уравнения переноса, из представлений (4.175)–(4.180) при кусочно-постоянной аппроксимации функции источника по z на отрезке [zk, zk+1 ] (4.165) получаем явную двухточечную схему:

Re+ ji = Re+ Dkj + Fkji A+ + Fk+1, ji Ckj ;

+ + Re+ Re+ (4.181) k+1, kji kj Im+ji = Im+ Dkj + Fkji A+ + Fk+1, ji Ckj ;

+ + Im+ Im+ (4.182) k+1, kji kj Re = Re ji Dkj + Fkji Ckj + Fk+1, ji A ;

Re Re (4.183) kji k+1, kj Im = Imji Dkj + Fkji Ckj + Fk+1, ji A Im Im (4.184) kji k+1, kj с положительными коэффициентами + 1 Dkj k + Ckj = A+ = + Dkj = exp (4.185), ;

+ kj 2tot, k j 1 Dkj k Ckj = A = Dkj = exp (4.186), ;

| | kj 2tot, k j k = tot,k zk.

Предложенный алгоритм построения конечно-разностного аналога сохра няет последовательность приближений и не перемешивает итерации: все расчетные формулы имеют вид рекуррентных соотношений, когда значения Re+, Re и Im+, Im на следующем уровне слоя на итерации с номером N вычисляются через значения перечисленных функций, взятых с предыдущего уровня слоя на той же итерации. В расчетные формулы параметр p входит в аргументах тригонометрических функций (sin и cos) – – это позволяет использовать любые значения ( px, py ) и получать ограниченные решения, но при этом Re и Im оказываются колеблющимися, и чем больше значения p, тем больше колебания. Поэтому в качестве критерия сходимости итераций предлагается проверять выполнение условия:

(N +1)+() N +() Xkji Xkji +() = max, (N +1)+() {k, j, i} Xkji 316 Глава 4. Комплексное уравнение переноса где – заданная относительная точность итераций, для слабо колеблющейся – функции 1/ Re+() 2 Im+() +() Xkji = kji + kji.



Pages:     | 1 |   ...   | 4 | 5 || 7 | 8 |   ...   | 15 |
 





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

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