ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2015, том 55, № 10, с. 1741-1755
УДК 519.634
НОВЫЙ ВАРИАНТ МЕТОДА ДИСКРЕТНЫХ ОРДИНАТ ДЛЯ РАСЧЕТА СОБСТВЕННОГО ИЗЛУЧЕНИЯ В ГОРИЗОНТАЛЬНО-ОДНОРОДНОЙ АТМОСФЕРЕ^
© 2015 г. Н. И. Игнатьев*, **, И. В. Мингалев***, А. В. Родин*, **, Е. А. Федотова***
(* 117997Москва, ул. Профсоюзная, 84/32, Ин-т космич. исследований РАН; ** 141700М.о., Долгопрудный, Институтский пер., 9, МФТИ; *** 184209Мурманская о., Апатиты, ул. Академгородок, 26а, Полярный геофизический ин-т Кольского научн. центра РАН) e-mail: mingalev_i@pgia.ru Поступила в редакцию 17.06.2014 г.
Переработанный вариант 22.10.2014 г.
Предлагается новый вариант метода дискретных ординат для расчета поля собственного излучения горизонтальной однородной атмосферы Земли и других планет. Особенность этого варианта заключается в том, что для решения возникающей в методе дискретных ординат системы линейных уравнений используется метод матричной прогонки. Этот метод является точным и максимально экономичным по объему вычислений, а также относительно несложным в программной реализации. Созданная авторами программа, в которой реализован предложенный метод, показала быстродействие примерно в 2 раза лучшее, чем быстродействие имеющегося в свободном доступе пакета программ DISORT. Библ. 10. Фиг. 3.
Ключевые слова: уравнение переноса теплового излучения, численный метод дискретных ординат, метод матричной прогонки.
DOI: 10.7868/S0044466915100117
1. ВВЕДЕНИЕ
Во многих физических приложениях возникает необходимость расчета поля собственного теплового излучения атмосферы Земли или атмосфер других планет и их спутников в приближении, когда атмосферу можно считать горизонтально-однородной, а собственное тепловое излучение зависит только от высоты и зенитного угла: например при расчете скоростей нагрева— охлаждения атмосферы за счет переноса собственного теплового излучения или при интерпретации дистанционных измерений различных атмосферных параметров. Поле собственного теплового излучения горизонтально-однородной атмосферы описывается одномерным по пространству уравнением переноса излучения с соответствующими граничными условиями. Для численного решения этого уравнения применяются различные варианты метода Монте-Карло и метода дискретных ординат, а также метод последовательных порядков рассеяния и метод сферических гармоник (см. [1]-[4]).
Метод Монте-Карло для рассматриваемой задачи успешно применялся, например, в [5], [6], а общие положения этого метода изложены, например, в [1], [2], [4]. Этот метод является приближенным. Он удобен для использования параллельных вычислений. Особенность этого метода заключается в том, что его точность пропорциональна 1/ JN, где N — число рассчитываемых траекторий фотонов. Кроме того, при наличии оптически толстых слоев, в которых альбедо однократного рассеяния очень близко к единице (например облачный слой в атмосфере Венеры), необходимое для достижения заданной точности число рассчитываемых траекторий фотонов может существенно увеличиться по сравнению со случаем отсутствия таких слоев. Метод последовательных порядков рассеяния также является приближенным. Он эффективен только в том
1) Работа выполнена при финансовой поддержке программы фундаментальных исследований Президиума РАН № 22, а также РФФИ (код проекта 13-01-00063) и гранта Минобрнауки № 11.G34.31.0074.
1741
случае, если в атмосфере отсутствуют оптически толстые слои (оптическая толщина более 0.5), в которых альбедо однократного рассеяния превышает 0.6—0.7. В противном случае итерации этого метода сходятся очень медленно.
Метод дискретных ординат для расчета собственного теплового излучения в узком интервале спектра в горизонтально-однородной атмосфере заключается в следующем. Вводится сетка по зенитным углам. Поле излучения разбивается на конечное число потоков, с каждым из которых связан фиксированный зенитный угол введенной сетки. Интеграл по углам, задающий источник рассеянного излучения, аппроксимируется линейной комбинацией потоков. Уравнение переноса излучения заменяется конечной системой обыкновенных дифференциальных уравнений, описывающих изменения с высотой интенсивностей излучения с заданными зенитными углами. Далее проводится дискретизация по высоте. Атмосфера предполагается состоящей из М горизонтальных слоев, в которых индикатрисса рассеяния и альбедо однократного рассеяния внутри каждого слоя берутся либо постоянными (но могут изменяться от слоя к слою), либо зависящими от оптической толщины (зависимость изменяется от слоя к слою). Осуществляется переход от системы дифференциальных уравнений к системе интегральных уравнений по высоте, связывающих интенсивности излучения в узлах сетки по зенитным углам на соседних слоях по высоте. После этого осуществляется переход к системе линейных алгебраических уравнений относительно интенсивности излучения в узлах сетки по зенитным углам и по высоте. Этот переход осуществляется с помощью аппроксимации интегралов по высоте в интегральных уравнениях аналитическими формулами.
Различные варианты метода дискретных ординат отличаются между собой способом аппроксимации интегралов по углам рассеяния и по высоте, а также методом решения системы линейных алгебраических уравнений. В пакете программ О^ОЯТ (см. [7]), который находится в свободном доступе, для аппроксимации интегралов по углам рассеяния используются квадратурные формулы Гаусса—Лежандра, а метод решения системы линейных алгебраических уравнений основан на нахождении собственных значений и векторов матрицы коэффициентов.
Также в свободном доступе имеется пакет программ 8НЭОМТ (см. [8]). Процедура 8НЭОМ представляет собой комбинированный метод сферических гармоник и дискретных ординат. В этой процедуре используется итерационный метод для расчета функции источника рассеянного излучения в узлах пространственной сетки, а угловая часть функции источников представлена в виде конечной суммы сферических гармоник. Этот процесс эквивалентен методу последовательных порядков рассеяния. Необходимое для достижения заданной точности число итераций возрастает при увеличении альбедо однократного рассеяния и оптической толщины.
В этой работе изложена новая модификация метода дискретных ординат для расчета собственного излучения горизонтально-однородной атмосферы, имеющая две особенности. Первая из них заключается в том, что расчетная сетка по зенитным углам может быть произвольной. Вторая заключается в том, что для решения возникающей в методе дискретных ординат системы линейных алгебраических уравнений используется метод матричной прогонки. Этот метод является точным и максимально использует структуру матрицы коэффициентов системы для уменьшения объема вычислений. Он является более экономичным и более простым в реализации, чем примененный в программе Э^ОЯТ метод решения, использующий вычисление собственных значений и векторов матрицы коэффициентов линейной системы, которая имеет большую размерность. В случае наличия в атмосфере слоев с сильным рассеянием и слабым поглощением (например слои облаков на Венере и на Земле) итерационные методы могут сходиться медленно и требовать выполнения большого числа итераций для достижения приемлемой точности решения. В этом случае предложенный в данной работе метод имеет преимущество в точности и скорости расчета.
2. ПОСТАНОВКА ЗАДАЧИ ДЛЯ РАСЧЕТА СОБСТВЕННОГО ИЗЛУЧЕНИЯ АТМОСФЕРЫ
Будем считать атмосферу плоской и горизонтально-однородной и рассмотрим собственное излучение атмосферы с частотой V, которое будем считать зависящим только от высоты над поверхностью и от угла между направлением импульса фотона и вертикальным направлением. Этот угол будем называть зенитным. Иногда его отсчитывают от направления вниз.
Введем обозначения: и — косинус зенитного угла, г — высота над поверхностью, zmяx — высота верхней границы столба атмосферы, в котором производится расчет поля излучения, Т(г) — температура атмосферного газа, а 1(г, и) — интенсивность излучения с частотой V и зенитным углом,
косинус которого равен u, на высоте г, Б(Т, V) = 2hv3c 2(ехр^у/^вТ)) — 1) 1 — функция Планка, в которой h — постоянная Планка, ^ — постоянная Больцмана, с — скорость света.
Уравнение переноса собственного излучения атмосферы в данном случае (см., например, [1]—[3]) можно записать в виде
л&и} = и) + (! - ю(г)) Щ Т(г))} V) + и), (1)
где _(г) и ю(г) — строго положительные, соответственно, коэффициент экстинкции (объемного ослабления) и альбедо однократного рассеяния атмосферного газа на высоте г для излучения с частотой V, Б(г, и) — перенормированная плотность источника рассеянного излучения с частотой V на высоте г и зенитным углом, имеющим косинус и. Эта плотность задается формулой
1 /2 л ^
Sz, u) = ^ i/(Z' Ц Jv(ф))
-1 0
dw, (2)
где w и u — косинусы зенитных углов до и после рассеяния, ф — разность между азимутальным углом
излучения до рассеяния и этим же углом после рассеяния, v(w, u, ф) = uw + cosфл/( 1 - u2)(1 - w2) — косинус угла рассеяния, a x(z, v) — индикатрисса рассеяния для излучения с частотой v на высоте z на угол, косинус которого равен v.
Оптической толщиной слоя атмосферы между верхней границей и высотой z называют взаимнооднозначно связанный с z параметр
zmax
T(z) = J CT(z)dz, dT = -CT(z)dz, T^max) = 0 , T(0) = Tmax.
z
Заменяя во всех функциях в уравнении (1) зависимость от z на зависимость от т, это уравнение можно записать в виде
иЩи) = I(T, u) - (1 - ш(т)) B ( Т(т), v) - S (т, u). (3)
dT
Уравнение (3) следует дополнить граничными условиями. Стандартным условием на верхней границе для собственного излучения атмосферы является равенство нулю направленного вниз излучения:
I(u) = 0, u < 0, т = 0. (4)
На нижней границе направленное вверх излучение складывается из рассеянного поверхностью падающего излучения и из теплового излучения поверхности, заданного выражением (1- Q(v))B(Tp, v) , где Tp — температУра пот&рхнюсти^ Q(v) — альбедо поверхности для изл:учения с частотой v. Если считать рассеяние поверхностью изотропным, то условие на нижней границе можно записать в виде
0
I(u) =
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.